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MAGNETIC BREAKDOWN IN A FINITE ONE -DIMENSIONAL MODEL 


by Gabriel Allen 
Levis Research Center 


SUMMARY 

Exact solutions have been obtained of the Schroedinger equation for a 
model approximating an electron in a one -dimensional periodic potential acted 
on by a transverse uniform magnetic field. In this model the periodic poten- 
tial is approximated by a finite chain (20 atoms long) of "periodic" square 
veils , and the parabolic potential due to the magnetic field is approximated 
by a parabolic square-veil potential. Procedures are presented for obtaining 
solutions for arbitrary combinations of the veil depth of the periodic poten- 
tial Vo, the magnetic field strength H, and the ratio of veil dimension (2v) 
to the atomic dimension (a), 2v/a. Eigenvalues and vave functions have been 
computed for the case 2v = a/2 for several values of Vq and H. The solu- 
tions obtained are valid for arbitrary field strengths, but in order to detect 
magnetic breakdovn effects in such a short chain, very large values had to be 
used. Although the values of H in the computations vere of the order of 
kilotesla (tens of megagauss), it has proved possible to associate these states 
vith their zero-field counterparts. The subroutines used in the computations 
are also included. 


INTRODUCTION 

The study of the Fermi surface of metals became greatly accelerated 
folloving the publication of the often quoted papers in references 1 and 2. 

The profusion of theoretical and experimental papers on Fermi surfaces, vhich 
began soon aftervards, continues to the present day. Consequently, there is 
nov available a veritable catalogue of quite detailed Fermi surfaces for a 
number of metals (ref. 3). 

The results of parallel theoretical calculations vere, for the most part, 
in such good agreement vith the experimental results that various second-order 
features of the structure (e.g., spin-orbit splitting) could be examined in a 
meaningful vay (refs. 4 to 6). One concept in vhich interest vas revived is 
magnetic breakdovn (refs. 6 to 10), a phenomenon in vhich the connectivity of 
a Fermi surface in a given direction may be changed in the presence of magnetic 
fields . 

In all of the treatments of the effects of magnetic fields on the electri- 
cal properties of solids, computations can normally be carried out only for the 


extreme cases of very small or very large magnetic fields. It was felt that it 
might prove enlightening to examine an exact solution of even a simplified 
model that embodied some of the important features of the interaction in solids 
of electrons with external magnetic fields. 

In this connection, in work done at the Lewis Research Center, exact solu- 
tions have been obtained for the case of an electron in a finite one- 
dimensional chain of "periodic" rectangular well potentials acted on by a 
uniform magnetic field in a direction perpendicular to that of the periodic po- 
tential. The periodic potential consisted of wells of width 2w separated by 
hills of width 2h = a - 2w, where a, the distance between atomic centers, is 
the period. The solutions are valid for arbitrary values of a, w, and h. 

Numerical values have been obtained for all of the bound-state eigenvalues 
and some of the eigenfunctions for a = 3 angstroms. Since the number of in- 
dividual computations was a monotonic increasing function of the length of the 
chain, it was necessary to keep this length down to 20 atoms. The magnetic 
fields were then chosen to be large enough so as to make the effective magnetic 
potential comparable to the 6 depth of the wells in the periodic potential. For 
a chain of this length (60 A), the magnetic fields were, therefore, of the 
order of kilotesla (tens of megagauss). 

It may be noted that the behavior of an electron subjected to such high 
fields in this model actually approximates that of electrons in real solids 
subjected to reasonable fields of the order of a few tesla. From a dynamical 
point of view, the field should be of such a size that the magnetic part of 
the potential becomes an appreciable fraction of the total potential over some 
part of the electrons path. In this model, the electron can travel a maximum 
distance of 60 angstroms before being scattered; therefore, fields of the order 
of kilotesla are required to build up the magnetic part of the potential in 
such a small distance. In real solids, on the other hand, the mean free path 
may be orders of magnitude greater than this, and much smaller magnetic fields 
could then accomplish the same objective. 

It has been proved possible, however, to associate each of these states 
with its zero-field state. Thus, by following the wave function for the zero- 
field state to its high-field state, something may be learned about the be- 
havior of the system for moderate or intermediate fields. 

The computations were performed on the IBM 7091 at the Lewis Research 
Center, and the Fortran TV subroutines used are described and listed in the 
appendixes . 


SYMBOLS 


A vector potential 

a distance between atomic centers, period 

constant defined by eq. (49) 
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c 

E 

Ejp 

E g 

e 

G 

g 


velocity of light 

parametrized quantity used like energy, see eq. (2l) 
Fermi energy 

zero-field energy gap in band structure 
charge of electron 



H 

— > 

H 


H n U) 

H z 

h 


h 


j 

kyj ^z 
m 


N 


n 

P 


o¥ 

b n' b n 


V(x) 


V M (x’ ) 

W x,) 


magnetic field strength 
magnetic field 

n^H degree Hermite polynomial of argument | 

magnitude of z component of magnetic field 

half width of hill 

Dirac h, Planck's constant/2n 

index specifying atom from center 

components of wave vector 

mass of electron 

integer defined by eq. (19) 

number of atoms in positive half of chain 

integer defining the number of energy level, e n 

momentum 

cyclotron radius 

determinants, see p. 38 

one -dimensional potential 

approximation to V^x'), see eq. (17) 

potential due to magnetic field, l/2 mo^x' 2 



Vp(x) ^VpCx 1 ) defined by eqs . (13) and (14), respectively 
Vq well depth of periodic potential 

Vq energy due to magnetic field defined by eq. (18) 

w half width of well 

x , y j z c oor dinat e s 


x* 

x 0 

r 

e 

€ n 

v 

*(?) 

CD C 

Subscripts : 
e 

i 

M 

mag 

max 

N 

P 

SB 

WB 

x,y,z 

a^|B 


x - x 0 

defined by eq. (6) 

defined by eqs. ( 30) and (3l) 

energy 

energy levels defined by eq. (9) 
integer defined by eq. (l9b) 
wave function 
cyclotron frequency 

even solution 

index specifying atom from center 

approximate magnetic 

magnetic 

maximum 

number of atomic distances from center of one end of chain 
counting central well as 0 

periodic 

strong breakdown 
weak breakdown 
coordinates 

defined by eqs. (CIS) to (C18) 
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Superscripts : 


h Mil 
w veil 

0 sign of appropriate g 
+ 

DERIVATION OF WAVE EQUATION 
Free Electrons in Magnetic Field 

The wave equation describing a free electron in a constant magnetic field 
may be written as 

^ (p + ! t(r) = ei(r) (l) 

■where it is the vector potential, p is the momentum, and e is the energy. 
The spin of the electron is neglected. 

— * — > * 

If the magnetic field H is given by H z k, where H z is constant and A 

is chosen in the Landau gauge (ref. 11 ), 

A=K z (0,x,0) ( 2 ) 

Then equation (l) may be written as 

2 i [ p x + (Py + | H z x ) 2 + P z j*U,y,z) = e*(x,y,z) (5) 

The substitution 

i|/(x,y,z) = A(x)expfi(yk y + zk z )l (4) 


will result in the following simplification of equation (3): 
-ft 2 d 2 A 


2m ^2 


s (*% + f ^ x ) 2 + s A - 


eTv 


Equation (5) is the equation of a harmonic oscillator centered about 


with frequency 


x 0 = 


cftM 

eH z 




e 

me 


H 


z 


(5) 

( 6 ) 

(7) 


5 



Thus, if x* = x - xo, equation (5) may be written in the more familiar form 

eV = 0 (8) 


-ft 2 d 2 A , 1 2 ,2, , 

-sr « + 7 mn> x' A + 

2m ^ 1 2 2 


From this form, it is known that 


ftt 2 k§ 

2 m 


e„ - 


ft 2 k| 


: n " 2m 


= cu c ft 


( n + 1) 


or (using eq. (7) ) 


= 


eE 

z z 


2 m + me ^ 


( n + i 1 ) 


(9) 


Thus, the allowed energy levels of a free electron in a constant magnetic 
field are given by equation (9). 


Thus, 


The term A n (x) will be a harmonic oscillator function in (eH z /cft)(x - x 0 ) . 


*n,k y ,k z = expjiCyky + (x - x 0 )]exp[^ ( x _ x 0 ) 2 ] (l0 ) 

where H n (|) signifies the n th degree Hermite polynomial of the argument 
A further discussion may be found in reference 12. 


Electron in One -Dimensional Periodic Potential 
and Constant Magnetic Field 

-> *5 an electron is acted on simultaneously by a constant magnetic field 

H = H z k and a one-dimensional potential in the x-direction, then equation (l) 
must be replaced by 




\2 i -» 

) + V(x) ^(r) 


= et(r) 


( 11 ) 


When the previous substitutions are used in equation (ll), it may be reduced to 


-ft 2 d 2 A 
2m qx2 


i L + eHz V 

5 x j 


2m 


ft 2 k 2 + V(x) 


A = eA 


( 12 ) 


The following assumptions should now be made: 

(l) The term V(x) is a periodic potential Vp(x) of period a such that 
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V P (x) = 0 
= V 0 


-w ^ x £ 
w £ x £ a 


•-} 


(13a) 


and 


Vp(x + na) = Vp(x) (13b) 

where n is any integer. Here w will be referred to as the well region. 

(2) The term xq is a whole number of atomic distances. Assumption (2) 
will have the effect of making Vp(x ! ) periodic with the same period a as 
Vp(x) since the center of a well region in x will coincide with the center of 
a well region in x 1 (equal to x - xq) so that 


V p (x T ) = 0 
= V 0 


-w £ X 1 £ w ^ 
w <; x 1 1 a - ¥ J 


and 


Vp(x l + na) = Vp(x* ) 

Then,, the equivalent of equation (8) can be written as 


^h 2 d 2 A 
2m 


dx 


A 1 2 .2, 

,2 + 2 * + 




A + Vp(x* )A = 0 


(14a) 


(14b) 


(15) 


The quantity l/2 mu^x' 2 , which will be called V ma g( x ')j may be considered to 
be a potential due to the magnetic field. Then equation (15) can be written as 


d 2 A 


dx 


,2 


22. L _ _ [y 

ft 2 r 2m L Vma S 


(x* ) + V P (x' ) >A = 0 




(16) 


Vmagfx') + Vp(x') 





Figure 1. - Form of potential due to constant magnetic field and 
periodic square well, w = a/4. 


The form of the potential Vmag( x 1 ) 

+ Vp(x') is shown in figure 1. 

In each interval, the potential is 
like a harmonic oscillator in that it is 
parabolic in x 1 . The general solution 
to such a potential problem is a sum of 
two Weber functions (ref. 13), each of 
which has special behavior at the origin 
and at infinity. In an actual harmonic 
oscillator, only one of these functions 
serves as a wave function because of the 
requirement that wave functions be well 
behaved at infinity. Since each interval 
in the present problem is finite. 
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however, neither of the Weber functions can be eliminated as an allowable solu- 
tion, and the usual linear combination of two independent solutions of the 
second-order differential equation must be used for the wave function. 

The form of the wave function is invariant throughout the entire range of 
x', the wave function for a given interval being distinguished entirely by the 
coefficients appropriate to that interval. 

It may be noted that the symmetry which seems to be physically apparent 
in this problem in the form of the potential is somewhat obscured by the mathe- 
matical form of the wave function. The Weber functions for -x are different 
from those for +x, and this fact must be taken into account when matching at 

x' = 0. 

There is, in general, no actual physical symmetry here. The term Vp(x' ) 
is symmetric about x = 0, whereas V m ag(x’ ) is symmetric about x = xo = 
-chky/eH z . The special assumption has been made that xq is an integral 
number of atomic distances. As was already said, in such a case the entire 
potential V^x') + Vp(x') is symmetric in x* . Clearly this assumption will 
only be satisfied for special combinations of ky and H z . 

However, the noncoincidence of the vertex of the parabola V ma g(x*) from 
the center of the well in Vp(x') (see fig. l) should not seriously affect the 
physical results as long as the resulting total potential does not differ 
greatly from the V ma g(x') + Vp(x'), which is the ordinate of figure 1. The 
depth of a well in the periodic part of the potential should be of the order of 
a few electron volts. The part of the potential due to the magnetic field 
turns out to be expressible as 

Vg< x '> “ 8-8X10-10 H 2(x - x 0 )8 jn- 


where H z is in tesla, and x - xo is in angstroms. 

The one -dimensional chain in the model is about 20 atoms long; therefore, 
magnetic fields large enough so that Vmag(x' ) contributes a few electron volts 
before the chain ends should be used. Since there are 10 atoms on each side of 
the center, this means V ma g(x , )« 8X10" 9 H§ electron volt at the ends of the 
chain so that H z « 1 kilotesla or 10 megagauss in order that V ma g(x') con- 
tribute 1 electron volt to the total potential. Because the potential varies 
as ( x f ) 2, V ma g(x 1 ) will be much smaller than 1 electron volt over most of the 
chain so that a deviation from symmetry should have a relatively small effect 
on the results. 


Approximation of Magnetic Potential by 
Parabolic Square -Well Potential 

Furthermore, a potential approximating VmagCx 1 ) should also not affect 
the results too drastically if it has the form 
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V M (x>) = m|< v o 

V M (x' + na) = V M (x') + (n 2 - v 2 )Vq 


where 


vi - -i- feY 

0 2m \ c / 


(17a) 

(17b) 

(18) 


H is in tesla, n is any integer, and m x i is an integer depending on x* so 
that 


m x i = v 


m x » = 0 O^x'^a-w 
va - w <; x 1 <L (v + l)a - w v any integer 


(l9a) 

(19b) 


The potential defined by equations (19) will be called a parabolic square- 
well potential. 


Define 


v(x') - ZB [v m (x') + Vp(x')] 




ft 2 k 2 ' 

2m 


( 20 ) 

( 21 ) 


where Vj^(x') is given by equations (l7a) and (l7b) and Vp(x' ) is still given 
by equations (l4). 

The fact that the total potential is not periodic means that the wave 
function and its derivative must be matched at each boundary between regions of 
constant potential. This means that the rank of the determinant used for the 

determinantal compatibility 
condition is at least twice 
as large as the number of 
regions of constant poten- 
tial in the positive half of 
n 2 Vq the system (taking advan- 
tage of symmetry). Thus, in 
order to keep the determi- 
nant down to a manageable 
size, the system must be 
kept finite. This is 
readily accomplished by add- 
ing the cutoff condition 



V[(N + l)a - w] = oo 


( 22 ) 
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where N is the number of atoms in the positive half of the chain. 

The form of V(x') can be seen by examining Vji(x r ) + Vp(x') which is 
shown in figure 2 (compare with fig. l). The substitution of Vm(x') for 
Vmag(x') has the advantage that the solutions in each hill or well region are 
now combinations of trigonometric or hyperbolic functions rather than Weber 
functions. The effect of Vm(x') is a kind of averaging of Vma g (x' ) in any 
region considered. The averaging could be improved by using (m2,+ m x , + i/3)Vq 

instead of m^tVQ. When considering the nature of the approximations already 
inherent in the model, however, the indicated improvement would be second order 
at best. 

Another factor that must be considered is the effect of the cutoff as 
given by the boundary condition in equation (22). It is clear from figures 1 
and 2 that the effect of the cutoff will be small for states lying well below 
Vo + N%o> since the wave functions would have small amplitudes when 
x' < (N + l) 2 a - w even if there were no cutoff. The energies for a few of 
the states investigated were not really low enough to avoid a forcing of the 
wave function to zero in the vicinity of x* = (N + l) 2 a - w. Most of the 
states examined are not greatly affected by the condition in equation (22), 
however . 


SOLUTION OF WAVE EQUATION 
Derivation of Wave Function 

When the quantities defined by equations (20) and (2l) are used, the wave 
equation for A(x' ) becomes 


+ [E - V(x')]A = 0 (23) 

dx'^ 

The solution of equation (23) is readily obtained. Denote the region from 
-a to a as the region. This is the region for which Vm(x') = 0. 
Furthermore, denote the interval by a subscript (starting with 0) and indicate 
by an h or w superscript whether the periodic part of the potential is 
"on" (hill) or "off" (well). Then in the O^h well, 

A^(x' ) = sin Y'qX 1 + D^ cos y^x r -w <1 x' ^ w (24) 

and in the O^ 1 hill, 

A^(x' ) = sin yj^x' + D^ cos yj^x' w <L x' <1 a - w (25) 

where the C's and D's are constants and 
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( 26 ) 


r w = 

'o 




(27) 


An analogous expression exists for the wave function in the 0^ hill for 
negative x ' ; that is, -(a - w) ^ x ' <L - w. 

In the n^ interval, 


A^(x' ) = CjJ sin Tnx' + cos y^x' na - w < x' < na + v 


>^(x' ) = sin yj^x' + cos y^x' na + w ^ x' ^ (n + l)a - w 


(28) 

(29) 


where 


and 



(30) 


(31) 


Equations (28) and (29) are of the proper form only for real positive 
y n f s. The other possibilities will be shown subsequently. 

Both even and odd solutions exist, in general. The even solutions will be 
discussed in detail and necessary modifications in procedure for odd solutions 
will be indicated where relevant. 

The form of the solution in the 0^ well becomes (subscript e denotes an 
even solution) 


A 0,e^ xT ) = D 0,e GOS r 0,e x1 -v ^ x l £ w (32) 

The solutions in the other intervals maintain the same forms as given by 
equations (30) and (3l), but pains must be taken to use the even coefficients 
and eigenvalues. Thus, 

fx’ ) = Q* s ± n yw x* + _ cos rif ^x T na - w < x 1 < na + w (33) 

n, e n,e n,e n,e n,e w 




n 



and 


A^x') = c i,e sin ^n, e x ’ + cos ^,e x ’ 

na + w^x 1 ^(n + l)a-w (34) 

where 

^n,e = e (35a) 

r%,e = - ||v 0 (35b) 

and 

E^a = E e - n 2 gvi (36) 

* & C 

The term E e is an eigenvalue of equation (23) belonging to an even eigen- 
function. 


It will be convenient to write the wave function in each interval in such 
a form that it is centered about one of the boundaries of the given interval. 
Dropping the e subscript for simplicity results in 







x* ) 

II 

cos y^x’ 


-V 

£ x' £ w 


(37) 

>£(*' 

) 

£ £! 
<C 

H 

sin 

[«■- 

(na 

- w) 

J + cos 

r w 

1 n 

[, 

“ (na 

- v) 













na - 

w £ x* £ na + w 

(38) 

a£(x’ 

) 

II 

> 
tl fT* 

sin 

n 

x’ - 

(na 

+ w)J 

+ b£ cos 
n 

r h 

* n 

x’ 

- (na 

+ w)J 











na + 

w <, X* 

< (n + l)a - w 

'(39) 


and 

An(x’ ) = AjJ sin Tw/ X ' - [(N + l)a - wjl (40) 


where the coefficients are all linear combinations of the C's and D's in 
equations (33) and (34), Bq = Dq e , the y's are given by equations (35), and 
N is the number of atomic distances from the center to one end of the chain 
counting the central well as 0. It should be noted that Aq is centered about 
x' =0; the arguments of and A^ are zero at the left boundary of the 

appropriate interval, and the argument of Ajj is zero at the boundary at which 
the potential becomes infinite and the wave function vanishes (which is the 
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right boundary of the last interval). 

As has been stated previously, solutions of the form shown are for real 
nonzero r ! s. If, for some eigenvalue e, x In some interval (say the j^) 
is zero, then in that interval 

Aj(x') = A j[ x '" (<i a ± w )] + B j (41) 

where the + or occurs according to whether the interval in question is a 
hill or well region, respectively. 


Finally, if X is imaginary in some internal, the trigonometric functions 
become the corresponding hyperbolic ones and the general form remains un- 
changed. Quantities g and G are defined as follows: 


s = r z 

(42) 

l r l 

(43) 


The form of the solution in a given interval will then be determined by the 
sign of g in the interval. (Note that gQ is always greater than 0.) Use 
will be made of the well-known properties sin ia = i sihh a and sinh ia = 
i sin a in the expressions that are used for A(x f ). Also, at this point the 
coefficients A will be redefined in such a way that the transition between 
positive and negative g values will be smooth. A summary of the different 
forms of A(x’ ) in various regions is 


Aq(x' ) = Bg cos Ggx ! 


where -w <£ x T C w, 

A? r 

— sinh x’ - ( 

G w L n 

n 


na - w 


)] + cosh G^[x' - (na - w)] 




A n[ x ’ " ( na - w ' 1 ] + B n 


— sin G^jjx' - (na - w)J + cos G^jx' - (na - w)J 


G w 
J n 


where 


na - w ^ x* ^ na + w. 


(44) 


g £<0 

4 = 0 

g£> 0 


(45) 
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A^(x ! ) = ^ 


^ sinh G^Jx' - (na + w)J + cosh G^jx' - (na + w)J 
A^^x' - (na + w)J + B^ 

sin G^jx' - (na + w)J + B^ cos G^Jx' - (na + w)J 


§£< 0 
gg = 0 


gh >0 




(46) 


where na + w < x 1 <(n + l)a - w, and 




j 5in t 

where Na + w ^ x' (N + l)a - w. 



S§ < 0 


• h = 0 


Sw 


g{j > 0 


(47) 


Eigenvalues. - The determinantal compatibility condition from which the 
eigenvalues may be determined can be obtained by matching the wave function and 
its first derivative at each of the boundaries. The details are given in 
appendix A. The determinant is large and unwieldy, but fortunately a system- 
atic procedure for evaluating it in general 
could be formulated. A subroutine was devised 
from which the results could be obtained on the 
IBM 7094. Figures 3(a) to (e) show the eigen- 
values for a chain of 10 atoms on each side of 
the Ota well plotted as a function of magnetic 
field strength (by using Vq). Each figure rep- 
resents a fixed periodic well depth (Vq). Note 
that the eigenvalues are actually e - ft2k§/2m 
rather than e (see eqs. (30) and (3l)). 


Wave functions . - The matching conditions 
that are used collectively to obtain the eigen- 
values can now be used individually in succes- 
sion to find the ratios of the coefficients in 
each interval to some given coefficient that 
will remain arbitrary except for normalization. 
The coefficient of the Qta well Bg has been 



(a) Even eigenvalues; V 0 , 1 electron volt. 
Figure 3. - Eigenvalues as function of Vq. 


14 



.20 


.16 


.12 


-o 

> 


.08 — 


.04 


0 


Breakdown 
boundary 
■ Blount 
- Strong 



(b) Even eigenvalues; Vq, 3 electron volts. 

Figure 3. - Continued. Eigenvalues as function of Vq. 


chosen as the arbitrary coefficient, and all of the other coefficients have 
been obtained in terms of it. The details are shown in appendix B. 


In order to compare the behavior of different states, it is convenient to 
use normalized wave functions. The normalization is shown in appendix C. 


MAGNETIC BREAKDOWN IN MODEL 

In order to see what magnetic breakdown means as applied to this model, 
the behavior of the system may be examined as the magnetic field increases, and 
an attempt to account for the changes in a reasonable way may be made. There 
are a few guidelines that may be laid out in advance without (it is hoped) 
prejudicing the interpretation. 

When there is no magnetic field, the system may be described as an elec- 
tron interacting with a periodic squar e-well potential in a box. There is a 
basic objection to the use of ordinary perturbation theory in estimating the 
effect of magnetic fields on the properties of laboratory sized samples 
(ref. 14). In fact, it is the periodic part of the potential that is commonly 
treated as a perturbation to describe the behavior of systems under the influ- 
ence of both a periodic potential and magnetic fields (refs. 7 to 10). Never- 
theless, the behavior of the system will change in a continuous manner as the 
magnetic field enters the picture. Speaking qualitatively, very small fields 
should have a relatively insignificant effect, and, as the fields increase, 
their effect should become more easily perceptible. Therefore, some measure 
may be sought that might be expected to indicate when the size of the field is 
such as to have readily discernible effects on the behavior of the system. 
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(c) Even eigenvalues; Vq, 5 electron volts. 

Figure 3. - Continued* Eigenvalues as function of Vq. 

Perhaps the simplest and most naive measure that comes to mind is a com- 
parison of the magnitude of the magnetic potential with V q, the depth of the 
well in the periodic potential. When equations (17) are used, this condition 
may he written as 

V M (x‘)=V 0 (48) 


The conditions described by equation (48), in which somewhere in the chain 
Vm(x') becomes as large as Yq, will be referred to as weak breakdown (WB). 

More explicitly, define 


G h = 



8.793984X10" 10 


ev 

(A) 2 (tesla) 2 


(49) 


and define n^ to he the value of v in equation (19) at weak breakdown. 
Then Vo is expressible as 


Yq = C H a 2 H 2 = i mo|a 2 


(50) 
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(52) 


Furthermore , 


n WB 


%B = 


C H H?a 2 


■1 

I ~ V 0 ~ 

V C H < ' :n WB a ^ 


(53) 


There is a second type of breakdown that may be called strong breakdown. 
In this situation, the change in Vm over a distance a equals Vq, so that 
the magnetic potential "washes out" the periodic potential. Thus, 


AV m (x') 

a sr- 


v 0 


(54) 


so that letting the subscript SB denote strong breakdown enables to be 

expressed as 


(%) SB = C H H 2 ( n SB a ) 2 


(55) 


n SB = 


V 0 

2C H H 2 a 2 


(56) 


and 


V ’ 

~ 2 ( 5 

2C H n SB a ^ 

Recently, a quite different criterion for breakdown has proved useful in 
explaining certain experimental results (refs. 7 to 10). This condition is 
sometimes called Blount breakdown and for an infinite crystal is commonly ex- 
pressed in the form 


E' 




&b c Ejp 


(58) 


where Eg is the zero-field energy gap in the band structure and Ep is the 
Fermi energy. 


e F> 


In this model, if it is assumed that an eigenvalue E plays the role of 
Blount's criterion becomes 


H 


8X10 3 ^2- 
E 


tesla 


(59) 
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or in terms of Vq, 


, 0.5 X E£ 

V 0 » — s. ev 




(60) 


It may be seen from figure 3 that actual gaps at zero field can be found, 
and these gaps may be used directly for Eg instead of resorting to first- 
order perturbation theory as is more commonly done (ref. 10 ). 


RESULTS AND DISCUSSION 


Description of Specific Model 

First, numerical values for the specific model used in the computations 
will be given. A value of 3 angstroms for a has been chosen as representa- 
tive of a large number of actual solids. It should be emphasized that the re- 
sults in the preceding section are valid for all values of a, w, h, and H z . 

In this report, the computations have, nevertheless, been limited to the case 

2w = 2h = a/2 = 1.50 A (6l) 

With this value of a , equation (18) is expressible as 


V 


0 = 


7.9X10" 9 h| — — - 

( tesla Y" 


Next the length of the chain was fixed by the time required for the sub- 
routines to go through a set of eigenvalues for fixed Vo and Vo. It turned 
out that a chain 10 atoms on each side of the 0^ well had a determinant of 
such a size that a set (with fixed Vq and Vq) could run in the maximum 
allowed time of 5 minutes. Thus, N was set equal to 10 in the computations. 


Eigenvalues 

The actual Fortran IV subroutines used in the computations are described 
in appendix D. The eigenvalues for various well depths are shown as a function 
of magnetic field strength (or Vq) in figure 3. For each Vq, computations 
have been carried out for a few values of Vq and these points connected by 
straight lines. (This procedure accounts for the kinks in the figures.) Fig- 
ures 3(a) to 3(d), however, show only even eigenvalues, while eigenvalues ob- 
tained from odd and from even solutions are shown in figure 3(e). These eigen- 
values are distinguished on the figures (where possible) by broken and solid 

lines, respectively. It may be noted that the chain in the model was long 

enough for the system to show a clear-cut band structure at zero field. As can 

be seen, the magnetic field shifts these zero-field levels by unequal amounts. 
Thus, the energy gaps in the band structure Eg (which are zero-field concepts 
in ordinary band -structure language) would be difficult to discern at large 
fields were it not for the fact that they were connected to the zero-field 
positions . 
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Some of the even and odd levels for higher fields are almost degenerate. 
These situations are denoted by a D at the energy in question. The actual 
computations that were performed using double -precis ion arithmetic show that in 
every case there was a nonzero separation between even and odd evergy levels. 
The pattern was always such that for a given Vq and Vq the lowest state is 
even, and then alternate odd and even states follow as far out as the computa- 
tions were carried. This is certainly to be expected and was only used as a 
rough check to see if any states were skipped in the eigenvalue search routine. 

Certain features to be expected if an actual band structure were observed 
were common to the eigenvalues of the system independent of the well depth and 
are shown for Vq = 10 in figure 3(e). First of all, there were always 21 
eigenstates in the first band. The lowest eigenstate was always even, and the 
eigenstate at the top of the first band was also even. The lowest state in the 
second band was therefore odd. 

Since it was intended to examine the three types of breakdown described 
in the preceding section, it was necessary to examine states in potentials for 
which strong breakdown had occurred. As a margin of safety for each value of 
Vq, computations were made up to eigenvalues 50-percent larger than the value 
of the pertinent (Vjq)gg. 

In this connection, it may be mentioned that, in order for the effects of 
strong breakdown to be manifest, the energy of the state being examined should 
be high enough so that the wave function has an appreciable amplitude in the 
region where the slope of Vj^(x*) is changing rapidly enough to wash out the 
effect of Vp(x') (see eq. (54)). Thus, strong breakdown will be said to occur 
in this model whenever the system is in an eigenstate such that E (Vjy[)g^. 

When equations (50) and (55) are used, the equation of the strong breakdown 
boundary is 


* - < 62 > 

0 

Both strong and Blount breakdown lines are indicated in the figures . The 
Blount boundary was plotted by using the actual gaps taken from the figures for 
Eg in equation (60). 

Figure 4 shows the shifting of the eigenvalues by increasing the well 
depth at a constant magnetic field. It may be noted that, for small well 
depths, the separation between successive eigenvalues does tend, as Vq goes 
to zero, to approach the constant value &d c (equal to 0.322 ev for the field 
chosen) as given by equations (7) and (9). It should be mentioned that the 
states on the E-axis which represent zero well depth are the actual values for 
a free electron in a magnetic field computed from equation (9). If the eigen- 
values for Vq = 0 in figure 2 are computed, the lower states (E S 7 ev) are 
rather close to those on the E-axis in figure 4, but the separation for the 
higher states becomes rather large (approximately 0.7 ev between the last two 
e igenvalue s shown ) . 
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Solution 


Even 

Odd 

Breakdown 

boundary 

Blount 
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0 Indicates even and 
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Figure 4. - Effect of well -depth at constant magnetic field of 2,8 kilotesla (28 megagauss). Values at 
Vq = 0 are eigenvalues for free electron in magnetic field obtained from equation (10). 


For the lower eigenstates, there seems to he a slight decrease in separa- 
tion as Vq increases . For the higher lying states in the first band, how- 
ever, the eigenvalues seem to cluster in degenerate odd-even pairs separated 
by about 1 electron volt. A comparison with figure 3 shows that the degener- 
acy is lifted after crossing the gap between the first and second bands. 

As the well depth increases, it is tempting to search for any tendency for 
the states to spread into a band (line broadening) as described in reference 
ence 10. However, this broadening arose from the degeneracy in the position 
of xq of each state at fixed Vq. By contrasty the computations in the 
present model were made at a fixed xq, so that there is no reason for the 
effect to show up in the model. 


Wave Functions 

The actual subroutines used in these computations are described in appen- 
dix E. In order to learn more about individual states, it is helpful to look 
at the wave functions. Wave functions have been computed and plotted for a 
large number of eigenvalues and some of the characteristics that were found are 
shown here. Advantage has been taken of symmetry so that the plots show only 
the positive half of the chain. 
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Most of the wave functions shown are for a rather deep well (Vq = 10 ev). 
Some wave functions for shallower wells have been computed and do not seem to 
differ qualitatively from these. Therefore, the discussion will be limited to 
this well depth unless the contrary is specified. 

The preceding discussion of the eigenvalues shown in figure 3 contended 
that the lines connecting the eigenstates could be interpreted as showing the 
change in the zero-field states as Vq increased. At a given Vq, the wave 
function for the lowest eigenvalue found (which is always even) has no nodes 
over the entire 20-atom chain; the wave function for the next one (which is 
odd) has one node at x* = 0, and the wave function for each successive eigen- 
value has one more node than its predecessor. The wave function corresponding 
to the i^* 1 eigenvalue thus has i - 1 nodes for any Vq. This fact makes it 
possible to follow the change in behavior of the ith state of the system in 
configuration space as Vq varies and permits a connection to be made between 
high-field and zero-field states. Of course, this test is not a sufficient 
one, since zero-field states could cross one another as the magnetic field in- 
creases. Thus, some additional factors will be considered in examining the 
behavior of the wave function to support the identification with indicated 
zero-field states in figures 3 and 4. 

In the discussion which follows, it should be noted that the number of 
nodes in the wave function for the 20-atom chain for even and odd solutions is, 
respectively, twice the number to the right of x* = 0 and twice the number 
plus 1. 

Another point requires some clarification. For a given eigenstate, a|J is 
determined by matching the wave function at the boundary between the last well 
and the last hill (where x* = 10a + w). Naturally, the vanishing of the de- 
terminant in appendix A is just the condition required to ensure that a|| 

would be exactly the same for this eigenstate if it had been computed by match- 
ing the dA(x')/dx' rather than A(x' ) itself at this boundary. If the deter- 
minant for a given eigenstate is sufficiently close to zero, then dA(x')/dx' 
is smooth everywhere. If it is not sufficiently close to zero, A(x' ) is smooth 
but )/dx' is discontinuous at x' = 10a + w and the amplitude of A(x' ) 

is inordinately large in this region. 

It turned out that it was not possible to match the boundary conditions 
with needed accuracy for all the eigenvalues . For this reason, a somewhat 
incomplete set of wave functions is. presented, in the sense that the same state 
can not always be followed for different magnetic fields . The figures show 
wave functions for states satisfying conditions as closely as were available 
under the circumstances. 

In this connection, it should be mentioned that the odd wave functions 
were most frequently inaccurate, and therefore, all of the figures show only 
even states except where otherwise indicated. As mentioned in the preceding 
section, the bottom of the second band is always an odd state. It is necessary 
to keep in mind that the lowest even state in the second band that is shown in 
several of the figures is actually the second state in the second band. 
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Location 

Middle of first band 
Top of first band 

Lowest (even) state in second band 
Middle of second band 
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Figure 5. - Wave functions for eigenstates at a magnetic field strength of 3. 1 kilotesla; V<J, 0.0765 
electron volt. 


The general features of the wave functions for the system will be examined 
in some detail. Figure 5 is typical and represents the system in a field of 
3.1 kilotesla (31 megagauss , Vq = 0.0765 ev). The lowest state in the figure 
is for E = 10.75 electron volts. This state is in the first band and has 
18 nodes. The amplitude of this state has a definite maximum at about 27 ang- 
stroms and is very small for x* < 22 angstroms. 

The next state shown is at the top of the first band (E = 12.34 ev). The 
general appearance of this wave function is quite similar to that of the pre- 
ceding state shown except that the maximum is even more pronounced and is 
shifted slightly to x T = 30 angstroms. 

The state following this one, although rather close to it in energy 
(E = 12.89 ev), demonstrates a sharply different character. For this state the 
amplitude is suddenly quite large in the vicinity of x 1 = 0 (although the 
maximums are not as pronounced as for the two preceding states ) and is very 
small for x 1 > 17 angstroms. 

The last state shown in figure 5 is for E = 15.29 electron volts and is 
in the middle of the second band. It shares with the preceding state a low 
amplitude near the end of the chain and a comparatively large one near the 
center . 


23 



Some semiquant itative statements may be made about this behavior instead 
of a full quantitative explanation. First of all, an explanation may be sought 
from the viewpoint that, in the first band, the system is largely controlled 
by the magnetic field with even the relatively deep (10 ev) "periodic M well 
having no more than a modulating effect on the basic free electron in a mag- 
netic field behavior. In this framework, it is proper to consider the cyclo- 
tron radius hk/rmo c as the primary parameter governing the motion of the elec- 
tron. For a field of 3.1 kilotesla (31 megagauss), this quantity is 

r « 6.35 -yfn + 1/2 A (63) 

Vhere n is the number of the state in equation (9). 

The state E = 10.75 electron volts is the 19th state, and so r should 
be about 27.3 angstroms, which is quite close to the peak of 27 angstroms for 
this state. The state at the top of the first band is the 23rd state, and, 
for it, r will be about 30.8 angstroms, which is still not far from the sharp 
peak at 30 angstroms . 

The next state, being at (or near) the bottom of the second band, may be 
expected to behave quite differently since it is on the other side of what 
would be the Brillouin zone for a truly periodic potential (see fig. 3(e)). It 
might be expected to behave more like a state near the bottom of the first 
band, so since it is the second state (the lowest state in the second band is 
odd), n may be set equal to 2 in obtaining an estimate of r. The resulting 
value of about 10 angstroms is not far from the actual region of large ampli- 
tude for this state. The final state in the middle of the second band behaves 
like a state that is less bound by the magnetic field than the others, and con- 
sidering the fact that its energy is comparatively high, this is not sur- 
prising. 

Another way of looking at the problem is to consider the motion in the 

k x - ky plane. States well below the top of the first band have orbits in 

k-space that do not come too close to the Brillouin zone. The state at the 

top of the first band has an orbit in k-space much of which is near the 

Brillouin zone boundary. On the other hand, the states near the bottom of the 
second zone have orbits in the second Brillouin zone and, in the reduced zone 
scheme, these are again not very close to the zone boundary. 

The same general type of behavior is shown in figure 6 for states at 
H = 6.2 kilotesla (62 megagauss). It is noted that r = 21.8 angstroms for 
n = 23 at this field and again this is very close to the position of the sharp 
maximum for the top of the band. The simple picture fails for the state in 
the middle of the band though, since no very clear maximum is present, and 
secondly, r would be about 19 angstroms, which is a rather low amplitude point 
here. Nevertheless, the overall behavior is quite similar to that shown in 
figure 5 . 

Some additional weight to this interpretation is furnished by figures 7 
to 9. Figures 7 and 8 show some odd wave functions that adequately satisfy the 
matching requirements. Figure 7 shows the states at H = 3.1 kilotesla 
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Figure 6. - Wave functions for eigenstates at a magnetic field strength of 6. 2 kflotesla; VA, 0.3061 
electron volt. 
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Figure 7. - Odd wave functions In first band; magnetic field strength, 3. 1 kllotesla. 
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Figure 8. - Odd wave function in first band; magnetic field strength, 4.4 kllotesla. 



Figure 9. - Even wave functions for well depth of 5 electron volts; magnetic field strength, 
5.9 kllotesla. 


(31 megagauss , compare with fig. 5). The wave functions would fit in properly 
with those shown in figure 5, and r would he 24 angstroms for the state in 
the middle of the band (not very good agreement) and 30 angstroms for the state 
near the top of the band (quite good agreement). 

Figure 8 shows an odd wave function in the middle of the first band for a 
field of 4.4 kilotesla (44 megagauss). The peak is at 21 angstroms, which is 
in very good agreement with r for this field strength. 

Some wave functions for a smaller well depth (Vq = 5 ev) are shown in fig- 
ure 9. Again the same general features as in figures 5 and 6 are exhibited. 

Thus, it would appear that the general behavior of the system described 
for figure 5 occurs under a variety of conditions . The picture used as a de- 
scription is too simple to be expected to apply uniformly to a more complete 
examination of wave functions for all combinations of well depth and magnetic 
field strength. Nevertheless, it appears to be somewhat useful in a limited 
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(d) States past Blount breakdown; V5, 0.306 electron volt; magnetic field strength, 6. 2 kllotesla. 

Figure 10. ** Blount breakdown In model; Vo, 5 electron volts. Dash-dot line, middle of first 
band; solid line, top of first band; dotted line, bottom of second band. 
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qualitative description. 

Blount breakdown, which was discussed in the section MAGNETIC BREAKDOWN IN 
MODEL can also be examined for its usefulness in predicting sharp changes of 
* behavior in the model. Figure 3 contains curves that show at which energies 
Blount breakdown is to be expected on the basis of equation (59) or (60). In 
examining Blount breakdown, changes must not really be sought in the behavior 
of states that are rather close to one another. Blount's criterion, as given 
by equation (60), actually states the breakdown occurs when E is "near" 
(2 Vq)1/2/e§; therefore, states near this quantity must be compared with others 
rather well separated from them. 

In order to test Blount breakdown in the model, a set of "Djj tested" 
eigenvalues (see appendix A) for a fixed value of H is required over a rather 
wide energy range. The results of the computations were such that such sets 
were available only for the well depth Vq = 5 electron volts. Thus, fig- 
ure 3(c) should be referred to in order to see where the states lie relative to 
the Blount line . 

It can be seen from figure 3(d) that all the states shown in figures 5 
to 9 are well below Blount breakdown. These figures, therefore, can furnish no 
information as to the validity of the Blount criterion for the model. 

Figure 10 shows wave functions for a well depth of 5 electron volts for 
states rather well before, near, and after Blount breakdown. Figure 10(a) 
shows states for Vq = 5 electron volts at H = 3.8 kilotesla (Vq = 0.1148 ev). 
These states share a feature with those states in figures 5 to 9 of all being 
below Blount breakdown (see fig. 3(c)). As in the aforementioned figures, the 
states flanking the energy gap demonstrate distinctly different behaviors. (It 
should be noted, however, that the roles of the states have been reversed in 
comparison with the earlier figures . The significance of this reversal is not 
clear at present). Next, states at H = 4.9 kilotesla (Vq = 0.191 ev) are 
shown in figure 10(b). From figure 3(c) it may be noted that the states flank- 
ing the gap, while still below Blount breakdown, are closer to it than the 
states in figure 10(a). The main feature of interest in this figure is that 
the difference between the flanking states is smaller than in figure 10(a). 

The amplitude of the state at the bottom of the second band is greater here and 
for A 1 = 1.5 angstroms is quite comparable to that of the state at the top of 
the first band. The latter state, in turn, has a greater amplitude near the 
end of the chain than the state in figure 10(a). The transition is completed 
at H = 5.8 kilotesla (Vq = 0.268 ev). As seen in figure 3(c), the Blount line 
passes through the gap at this field, so both flanking states are actually near 
Blount breakdown. These states are shown in figure 10(c). It may be noted 
there that the differences between the flanking states are greatly diminished 
in comparison with cases shown below the Blount line. Incidentally, the ampli- 
tudes of the flanking states are now in the same relative position as those at 
Vq = 10 electron volts, which were all below Blount breakdown. Finally, fig- 
ure 10(d) shows states past Blount breakdown. 

Thus, it would appear that the periodic part of the potential exerts a 
strong influence on the system for states well below Blount breakdown, the 
states corresponding to the ones flanking the first gap in the zero-field band 
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structure consistently exhibiting significantly different behavior. On the 
other hand, as Blount *s criterion begins to be satisfied, the periodic poten- 
tial has a much weaker effect, and the changes between states flanking the 
energy gap become far less pronounced. 

The remaining type of breakdown, strong breakdown, has not been investi- 
gated in detail for reasons that would seem to be obvious. The strong break- 
down line in figure 3, for the most part, involves states that lie so far above 
the bottom of the second band that they would lie in the continuum were it not 
for the box property of the chain. A few wave functions in this region were 
computed and showed only the typical behavior characteristic of "free particles 
in a box." These wave functions, therefore, would not appear to contribute 
markedly to an understanding of breakdown in the model. 


SUMMARY OF RESULTS 

A one -dimensional model has been examined in an attempt to follow magnetic 
breakdown in some detail. The model is that of an electron in a uniform mag- 
netic field H z and a one -dimensional chain 20 atoms long of periodic square 
wells with infinite potential at each end. Exact solutions have been obtained 
for both the eigenvalues and the wave functions for arbitrary values of well 
depth, well width, atomic separation, and magnetic field strength. Computa- 
tions have been carried out for several well depths from 1 to 10 electron volts 
and for several magnetic field strengths. Although the magnetic fields were of 
the order of tens of megagauss, it has proved possible to associate these 
states with corresponding zero-field states. 

Three types of breakdown were considered for applicability to the model. 
The simplest type was one in which the magnitude of the magnetic part of the 
potential was equal to the depth of the "periodic" part of the potential. 
Another type was one in which the field was so large that the increase in the 
magnetic part of the potential over an atomic distance was equal to the well 
depth. The third type of breakdown examined was Blount breakdown in which 
ho) c E^/E| is compared with 1. 

An examination of the wave functions for various conditions showed that 
below Blount breakdown the system is controlled by the periodic part of the 
potential and the accompanying band structure. The most persistently apparent 
characteristic of these results is that at a fixed magnetic field there is a 
sharp change in behavior of the wave functions in going from a state corre- 
sponding to the top of the first band in zero magnetic field to a state corre- 
sponding to the bottom of the second band in the zero-field situation. This 
behavior seems to admit an interpretation in terms of a change of orbits in 
k-space between states on either side of a Brillouin zone. 

When the magnetic fields become large enough to cause the system to under- 
go Blount breakdown, the aforementioned differences between states flanking the 
zero-field energy gap become attenuated to a large extent. Consequently, in 
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this case, the periodic part of the potential plays a far less decisive role in 
determining the behavior of the system. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, October 1, 1964 
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APPENDIX A 


DERIVATION OF EIGENVALUES 


As mentioned in the section SOLUTION OF WAVE EQUATION, Bq is chosen as 
arbitrary. Then the requirements of continuity of A(x') and (d/dx')A(x') are 
satisfied by matching the wave function at every boundary. The boundary be- 
tween the Oth hill and the 0 th well occurs at x' = w. Thus, 



x’^w 



(Al) 


and 




(A2) 


Substituting equations (44) to (47) into equations (Al) and (A2) results in 


Bq cos wGq = Bg 


(A3) 


and 


- B 0°0 sin wGq = A§ (A4) 

Both equations (A3) and (A4) are independent of the value of gQ. The remain- 
der of the matching equations are considered next. At the boundary between the 
n^ well and the n^ h hill (this boundary is the right boundary of the well 
region and the left boundary of the hill region), x 1 = na + w and it is there- 
fore required that 


A^(na + w) = A^(na + w) (A5) 

and 

-AL- A^(na + w) = ~r A^(na + w) (A6) 

By examining the form of A^(x’ ) and A^(x' ) from equations (45) and (46), 
it is noted that A£(x' ) involves functions of x' - (na - w), whereas A£(x' ) 
is expressed as a function of x 1 - (na + w). Thus, A^(na + w) will be a 
function of 2w and A^(na + w) will be a function of 0. Therefore, the match- 
ing equations will take the form 


% 
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and 


(A^/G^)sinh 2wG^ + Bj[ cosh 2wG^ 
A^(2w) + B£ 

(A^/G^)sin 2wG^ + B^ cos 2wG^ 


cosh 2wG™ + B^G™ sinh 2wG^ 


A v 

n 


A^ cos 2wG^ - B^G™ sin 2wG™ 


,W tdVa¥ 


^W 


3-W < o 

=>n 


jw 


■\ 


Sn = 0 ^=B* 


p* w > 0 
& n 






.w < o 




= 0 
hi 


> = A* 1 
n 


Sn> 0 




Both equations (A7) and (A8) are independent of the value of g^. 


(A7) 


(AS) 


In a very similar way; the form of the matching equations at the boundary 
between the n™ well and the n^* 1 + 1 well (where x 1 = (n + l)a - w) can be 
readily expressed as 


(A^/G^)sinh 2hG^ + B^ cosh 2hG^ < 

A h (2h) + B h g h = 

n K 7 n & n 

(A|j/G^)sin 2hGj^ + B^ cos 2hG^ g^ > 

A^ cosh 2hG^ + B^G^ sinh 2hG^ g^ < 

A h rrh — 

n ®n 

A^ cos 2hG^ - b£g£ sin 2hG^ g£ > 


"\ 


> = B w _ 
' n+1 


> = A n + 1 


J 


(A9) 


(A10) 


Both equations (A9) and (A10) are independent of the value of g^ +1 - (The 
simple form of the right side of these equations is one of the main reasons for 
centering the wave functions as in eqs. (l5) and (16)). 

These equations represent all of the conditions except the matching at the 
last boundary in the chain. In the (N + l)^ hill; use must be made of equa- 
tion (47) for (x 1 ) at x* = Na + w. This value of x 1 will make AjJ(x ! ) a 
function of -2h. Therefore; matching at this last boundary yields the follow- 
ing: 
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• sinh 2hc£ 
P h K 


-7 sinh 2wG n + B H COsh 2wG N = V A N (2h) 


4 h 

- 3 sin 2hG* 
V, °N 


r 


Am 


.sinh 2hG^ 
,h N 


Ato- 

— sin 2wG^ + BjJ cos 2vGjJ = ^ 

% 


aJJ 

- - sin 2hG* 


aJ 

- 3- sinh 2hGjJ 

% 


- A^(2h) 

A? 

- — sin 2hd 


TJ 


4 <0 


«N 


= 0 


K<° 


4 >0 

J 


< 0 


Ag(2v) + B$ = J - Ag(2h) g§ = oUjJ = 0 


<4>° 


6w < 0 


Sw = °f g ?5> 0 


g W > 0 

J 


A*J cosh 2hGj| 

cosh 2wG^ + bJJgJJ sinh 2wGjJ = ^ A^ 

cos 2hG N 


cosh 2hGjJ 


A N " i 

<4 cos 2hGjJ 


cosh 2hG^ 


AjJ cos 2^0^ - B^G^J sin 2wG^J = ^ 


v.a£ cos 2hG N 


s|j < 0 

g N = 0 ^ g N < 0 

4 >0 - 


h ^ ^ 
gjf < 0 


g N " 0 ^ g K “ 0 


g N > 


4 < ° A 


g w = 0 


r g S > 0 
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In order to obtain the determinantal compatibility conditions from these 
equations , it will be convenient to introduce some notation that will enable 
the three possible expressions for the matched wave function (which depend on 
the values of g) to be condensed into a single expression. 


A C preceding an A or B will denote the coefficient of that A or 
B in the appropriate matched wave function A, and C'A or C'B will denote 
the coefficient of the A or B in the appropriate matched dA/jx'. A sub- 
script L or R on the C will denote whether the matching is at the left 
or right boundary, respectively, of the region in which the wave function is 
operating. Thus, w is the right boundary of the 0 th well and the left 
boundary of the hill, so that equations (A3) and (A4) can be written 


b o c r b o - B S - 0 
b o°r b o - A o - 0 


where 


„ ^w «w 

C'rBq = cos wG 0 


C R®0 = - G o sin wG o 


gY always > 0 


(A13) 

(A14) 

(A15) 

(A16) 


In a similar way, it is seen from equations (A7) and (A8) that it is possible 
to write 


+ b M - - o 


A^C^A^ + B w c'B w - A h = 0 
n K n n R n n 


(A17) 
(A18 ) 


so that 


r 


+ aW 


CIA 


R n 


sinh 2wG£/G£ 
2w 

sin 2wG^/G^ 


(A19) 


pS-oW 

C R B n 


0 

G R + ^ = < 


r cosh 2wG^ 


V. 


cos 2wG^ 
n 


(A20) 
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I 


0 

c r\ = < 


~Cr* sinh 2wG w 
n n 


-G^ sin 2 wG£ 


The superscripts + on the C's refer to the sign of the appropriate 
Equations (A9) and (A10) may be written in the same way. Thus, 

+ b M - B 2 + i - 0 

+ - A n + 1 * 0 

where 


r 


o 

C+A h = J 
R n ^ 


sinh 2hG^/ G^ 


2h 


sin 2hG^ 1 /G^ 1 
n' n 


V. 


''cosh 2hG^ 


C+B h = C ' A h = < 
R n R n ^ 


cos 2hG^ 
n 




^G 11 sinh 2hG h 
n n 


c e + b S - < 


-G^ sin 2hG h 
n n 


Finally, equations (All) and (A12) become 


(A2l) 

(A22) 

(A23) 

(A24) 


(A25) 


(A26) 
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(A27) 


A I°R A I + ■ A N C ’L A N - 0 

+ B S C R^ " A R C L A N = 0 (A28) 

where CpA^J, C R B^, C R A W^ and C R^J have the same form as equations (A19) 
and (A20) with n = N and where 


-sinh 2hG^/GjJ 




-2h 


-sin ShGjj/Gjj 


0 

C' + A h = ■( 
L N ^ 


cosh 2hGn 


N 


cos 2hG: 


h 


N 


(A29) 


(A30) 


The usual argument is now invoked, which says that if equations (A13), 
(All), (A17 ), (A18 ) , (A22 ), (A23), (A27), and (A28) are to hold simultaneously, 
then the determinant of all the coefficients must. vanish. This is the determi- 
nantal compatibility condition that is being sought. Equation (A3l) shows this 
determ inan t, which will be denoted by %. The columns are labeled according 
to the A or B whose coefficient appears in %. It will prove convenient 
to start with the equation for the last boundary and work toward the 0 th well. 






a£l 

a N-1 


b n- 


+c l a w - c M -W 0 
+c l a n ~ c A " c r*n . 0 


0 0 
0 1 
0 0 
0 0 

0 0 
0 0 
0 0 


-vS-i ~ c A-i 

- C X-1 - C E^-1 


BJS-l 
0 
0 
0 
0 

‘Vk-i - c r b k-i 

- C A-1 

0 0 

0 0 

0 0 


A S Bg 


-¥o -°r b 0 

Vo -^=0 


0 

0 

0 

0 

0 

0 

0 

0 

-cos vGY 


1 0 G v sin wG w 

0 0 


(A31) 
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The determinant is(4^+2)x(4N+2)« It will be evaluated in the 

usual way be expanding it in minors. By referring to (A3l), it is seen that 


~ C A oooo 

0 1 " c r a r-i " c r^-i 0 0 

1 0 -C^.! -C^B^ 0 0 


0 


% ” c l a tJ 


0 0 
0 1 


1 _C R A N-1 " C R e ®-1 ' 

0 “ C R A 5- 1 ^R^-l • 


0 0 
0 0 
0 0 
0 0 
0 0 


0 

0 

0 

0 

0 


0 0 0 
0 0 0 
0 0 0 
0 0 0 


0 0 0 
0 0 0 
0 0 0 
0 0 0 




-C A n 

°r a o 

0 

1 


- c r b o 


0 


- c r b S 0 

X -cos wGq 


0 Gg sin wGg 


C L A ^ 


- g r^ - c r^ 


“ c r a n " c r^5-i 


-^-1 




0 0 

1 " C R A N-1 ‘ C R B W-1 

° - C R^-1 - C R B N-1 

0 0 0 


0 0 
0 0 
0 0 
0 0 
0 0 

- c r a $ - c r b £ 

“ c r a o - c r b o 
0 1 


0 

0 

0 

0 


(A32) 


-cos wG 


,w 


0 G” sin wGq 


It may be observed that the determinants which are the coefficients of 
and C^Al}, respectively, are identical except for their first rows. The coef- 
ficient of c l A W is denoted by Sjp while the coefficient of C^-A^ is denoted 
by s3. An attempt at descriptive notation is being made here since and 

are ; respectively, determinants in which the upper right terms are “CrA^ 
and -C R A^. With this notation 

% = = 0 (A33 ) 
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Nov let Sjj be reduced to see if a pattern can be found vhereby the 
entire determinant can be evaluated readily. Expanding Sjj by minors yields 


1 

o 

o 

0 



6 


- c RBW-1 0 0 

"^^-1 0 0 

0 1 - c R A 5- 1 - C R^-X 

1 ° - C R^-1 - C R^1 

6 6 6 6 


0 0 0 


o 


o 


ooo 


0 0 


0 0 0 0 


0 


0 0 
0 0 
0 0 
0 0 

■ c r a o ‘°r b o 
■°r a o ~ c r b o 
0 1 
1 0 


0 

0 

0 

0 

6 

0 

-cos wGq 
Gq sin wGq 


-cX 0 0 0 0 

1 - C R A W-1 - C R^-1 0 0 

00 1 “^N-l “ C R^[-1 

° 1 ° “V^-l - C X-1 


0 0 
0 0 
0 0 
0 0 


0 0 0 
0 0 0 
0 0 0 
0 0 0 


0 0 
0 0 
0 0 
0 0 

Vo -°r b o 

'Vo "Vo 
0 1 

1 0 


0 

0 

0 

0 

0 

0 

V 

-cos vGq 
Gq sin wGq 


(A34) 


Expanding each remaining determinant by minors again yields 

Sr = " C R®N S N-1 (A35) 

where and S^_j_ have a readily deducible connotation analogous to sJJ 

and SjJ, respectively. It should be noted that the determinants 

and have a structure much like that of SjJ and SjJ; they merely begin 

with " c r a r_]_ and - “ C rAr_i instead of - c rAJn and ‘ -C R A 1F respectively, and 
are of dimension two less than SjJ and SjJ. 

Similarly, expanding SjJ by minors twice yields 
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<S = “V* 


1 

0 

0 

0 

0 

0 

0 

0 


-Vh- 1 "VS-l 0 


- C R^-1 


0 

1 

0 

0 

0 

0 


-C B 


R N-l 
1 


0 

0 

0 

0 


0 

0 

0 

0 


- C R^-1 - C R^-1 

- C R^-! - C R^-1 


0 

0 

0 

0 


0 

1 


- c r a “ - c r b S 
“ C R A 0 "°r b o 


1 

0 


0 

0 

0 

0 

0 

0 


-cos 

G-Y sin wG^ 



- C R^ 0 

0 0 0 

, . . 0 

0 

0 





1 - C A - 1 

-°A-i 0 0 

, . . 0 

0 

0 





0 0 

1 - C R^-1 - C R B S-1 • 

. . . 0 

0 

0 




+ 

0 1 

° Viu - c r^.i ■ 

. . * 0 

0 

0 



(A36) 


0 0 

0 0 0 

■ • • - c r a o 

- c r b § 

0 





0 0 

0 0 0 

• • • - c A 

- C R B S 

0 





0 0 

0 0 0 

. . . 0 

1 

-cos wGq 





0 0 

0 0 0 

. . . 1 

0 

Go sin vGq 




Thus, 


q¥ _ p aVo h 

b N ~ TrVlT-1 “ 

p ■rW oh 





(A37) 

Since the form of 

SjJ. 2 . and is 

so similar to 

the form 

of 

oW 


and S£, respectively., 

it is clear that a general procedure now exists 

for 

expanding by minors in successive steps. Thus, 








s l - -VS6.1 - 

n -rWq 11 
°R B n s n-l 





(A38) 



SW = -C'A^S 11 . - 
n R n n-l 

C lM-l 





(A39 ) 

and 


s£ " - 

C E^ 





(A40) 
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- -c^55 - 


(A4l) 


The process continues in this fashion until Sq and Sq are reached. 
These have the special form indicated by the three columns and four rows in the 
lower right corner of all of the previous determinants. Thus,, 


or 


and 



-C A* 1 

TO 

0 

.1 


- c r b S 0 

1 -COS wGq 

0 Gq sin wGj 


Sq — -Cj^AqGq sin wGq + C-^Bq cos wGq 



-^R^O " C R B 0 0 

0 1 -cos wGq 

1 0 G^ sin wG^ 

0 0 


(A42) 


so that 

S§ = -C^aJgJ sin wGq + cos wGj (A43) 

Actually, equations (A42) and (A43) can be made consistent with equa- 
tions (A38) to (A4l) if the following identification is made: 

Sq = -cos wGq (A44) 


s£ = Gj sin wGq (A45) 

In going over what has been done, it may be noted that a procedure has 
been developed for evaluating D N by starting at the lower right corner and 
working toward the upper left corner. Since this is a somewhat unusual proce- 
dure, there may be some merit in summarizing it. 

The terms Vq and Vq are fixed and then the eigenvalue for such a 
V(x' ) is determined as follows: A trial value of E is chosen, Gq is 

evaluated, and then Sq and Sq are obtained from equations (A44) and (A45 ) . 
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These values are then substituted into equations (A40) and (A4l) and by succes- 
sive application of equations (A38) to (A4l) all of the S’s and S's through 
SjJ and sg are obtained. Then the quantity in equation (A33) is found, 

Djj = C L A]pjj - C l AjJs]J, and a decision is made as to whether. it is sufficiently 

close to zero. If it is, the trial value of E is chosen to be an eigenvalue 
for the fixed Vq and Vq assumed. If the 3% so formed is not close 
enough, the entire procedure is repeated with a new trial value of E until 
a sufficiently small is obtained. Actually, there are many eigenvalues 

for each given Vq and Vq, as can be seen from figures 3 and 4. 

It is readily seen that for odd A(x' ) the only change that need be made 
in the aforementioned procedure is that Xq(x’ ) = Ag sin Ggx’/Gg and that Sg 
and Sq will be given by 

Sq = -sin wGq/Gq (A46) 


Sq = -cos WGq 


(A47) 



APPENDIX B 


DERIVATION OF WAVE FUNCTION 

Once having obtained an eigenvalue, the unnormalized wave function for 
that eigenvalue is readily obtained from the matching equations shown in 
appendix A. Since Bq is going to be used as the arbitrary amplitude, there 
immediately results (from eqs . (A13) and (A14)) 

A o - b o c b b o - <- G 0 sln 

b{$ - BqCrBq = (cos vGq )Bq 

Thus, 

A S/ B S ■ 3ln < 

Bq/Bq = cos wGq 

Equations (Bl) and (B2) are then substituted into equations (A22) 
and (A23) with n = 0, the results substituted into equations (A17 ) and (A18), 

and these equations used successively until AjJ/Bq and BjJ/Bq are obtained. 

At this point A^J can be obtained from equations (A27) or (A30). If the value 
of Djj used to find the eigenvalue was sufficiently small, then the values of 

A^/Bq obtained from these two eq.ua tions will be close enough. If % were 
actually 0, then the values of a|j using the two equations would be identical. 

The set of coefficients of A(x') in each interval having been obtained, 
the (unnormalized) wave function at any given value of x 1 can be computed. 

The normalization of A(x') would allow a comparison to be made between the 
relative probability densities of eigenfunctions belonging to different eigen- 
values in any region. This procedure is carried out in appendix C. 

For odd A(x'), An will be used as the arbitrary amplitude, and instead 
of equations ( Bl ) and ( B2 ) , there will be 


a q/ a o = cos(vG’") 

(B3) 

B B /A» = sin(«Gp/0“ 

(B4) 


The remainder of the procedure is unchanged. 


(Bl) 

(B2) 
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APPENDIX C 


NORMALIZATION OF WAVE FUNCTION 

The wave functions will be normalized in the usual way; that is, the 
otherwise arbitrary constant Bq is determined so that 

| A(x’ ) [ 2 dx* = 1 (Cl) 


As has been seen, A(x’ ) assumes a different form in each interval and, 
within a given interval, for each sign of the appropriate g for that inter- 
val. Whether the solution be even or odd, it will still be true that 




Furthermore, the right side may be broken up in the following way (taking 
account of the fact that A(x ! ) = 0, x* ;> (N + l)a - w): 



The A(x' ) in each interval may be replaced by its special form as given in 
equations (45) and (46), and the normalization condition obtained in the form 



In appendix B, a procedure was displayed for obtaining each of the coef- 
ficients A^, B^, A^, and in terms of Bq. It is clear then that, if all 

coefficients be expressed in terms of B^, the wave function in each interval 
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can be written as B^f^x' ) or Bgf^(x')> where fg(x' ) and fg(x' ) are the 
forms of T'g(x') and Ag(x'), respectively, after the Ag, Eg, Ag, and Bg 
have been so expressed. 


If 


*z=r i f j( x, )i 2 <3x ' 

•'O ~ 


/*( n+l)a-w 

*g = / |fg(x ')| 2 ax' 

•nia+w 


(C5) 


(C6) 


and 


■jji w _ 

F n ~ 


/•na+w 


/ |fg(x')| 2 dx’ 

there results at once (using eq. (C4)) 

1/2 


-dW | 2 _ 

B 0l " 


F 0 + 


N N 


£ *£ + £ F n 
n=0 n-1 


(C7 ) 


(C8 ) 


The remaining task is to find the F's. The form of fg(x' ) is always 


cos G^x 1 ; therefore. 


rjiW 


✓»w 


cos 2 Gqx' dx' = | | w + 


sin 2vG ( 


.V 


0 


2Gq 


(C9) 


Each of the remaining F*s may have several forms depending on the value of 
g that goes with the f in question. It is readily ascertained that the 
change of variables 


v = x* - (na - w) for a veil region (G10) 

ja=x* - (na + w) for a hill region (Cll) 

will simplify the algebra involved in evaluating the F T s. With these changes, 
there remains only the evaluation of the following integrals: 
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3 ^ 


(C12a) 



(C12b) 


(C12c ) 

(C13a) 
(C13b ) 
(C13c) 

(C14a ) 
(C14b) 


(C14c) 



Equations (C12) to (C14) are the forms involved when g^ is negative , 
zero, and positive,, respectively. Corresponding well integrals are exactly the 
same with w replacing h each time the latter occurs. 


If the a 1 s and p’s 
Bq so that 


are defined as the ratios of the A's and B’s to 


W-nW _ A w 

a n B 0 ~ A n 


p w B v = B v 
K n 0 n 

-h-DW _ «h 
a n B 0 = ^ 


p h B^ = B h 
n 0 n 


(C15) 

(C16) 

(C17) 

(CIS) 


then 

^ = |ag| 2 [a] + |ph| 2 [b] + («£*?£ + =gpf)[c] (C19) 

where [a], [b] , and [c] are obtained from (C12), (C13), or (C14) depending on 

sS- 


Similarly,, 

Fg = |ag| 2 U] + | PK| E Cb] + («£*(£ + ogpg*)[c] (C20) 

where [ a], [b], and [ c ] are obtained from the equations corresponding to (C12), 
( C13 ) , and (C14), respectively, for a well region. 

By choosing the forms in equations (45) and (46) as has been done, all of 
the A's and B's will be real (and thus also the a‘s and (3*s). Then 
equations (C19) and (C20) may be written 

1* - ( a n) 2 [a] + Og) 2 [b] + 2^tc] (021) 

F~ . (a”) Z (a] + (^) 2 [bj + 3a^3“[o] (C22) 

Wow and can be written in full: 
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This completes the normalization. All of the quantities needed in equa- 
tions (C8) are given by equations (C2l), (C22), and ( C9 ) . If they be substi- 
tuted in the right side of (C8) and the square root of the result taken, then 
the value of BS so obtained will cause equation (Cl) to be satisfied. 

It should be noted that the results hold both for even and odd eigen- 
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functions, the only differences between these cases being that, first, a 
ficient other than B$ will have to be chosen as arbitrary (since Xq(x 
be of the form sin G^x'/G^j), which will make -Fg take the form 


F w _ 
0 


w 


(°o) 


1 - 


sin 2 wGq 


2wG, 


,w 


Secondly, the resulting a's and P's will be different for odd A(x’), 


coef- 
) will 

(C25) 
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APPENDIX D 


FORTRAN IV PROGRAM FOR COMPUTATIONS OF EIGENVALUES 


General Description 

This program written in Fortran IV language (actually a converted For- 
tran II program) was used to compute the eigenvalues shown in figures 3 and 4. 
It is arranged to operate -with the Levis Research Center 7094 monitor system. 
It consists of a main routine and four subroutines. The main routine and the 
first subroutine are also used in the program to compute the wave functions. 


Before listing the individual routines , it "will be useful to note the 
following: 


By using equation (56), ng-g may be written 


n 


SB 




2V 


0 


(El) 


so that 


(%) 


m'sb 


(-> TT I 

= n SB V 0 = 


V 0 

(2Vq) 2 


(D2) 


As indicated in the section RESULTS AND DISCUSSION, the computations have 
been carried out to values 50-percent greater than (V^)g^. Thus, the computa- 
tions are carried out to a value of H large enough to make V M (x')~ 3(V m ) sb /2 
or 


3 

2 





n2v 0 


from which a maximum value of 


Vq was chosen. 


Thus, 


(O 

v max 


v c 

2N 


(D3) 


(D4) 


Normally, this quantity was computed by the first subroutine VOPFIX for 
each value of Vq, and then computations were made for various multiples of 

(l/8)"kk of this value. For some purposes, computations were desired for a 
fixed (Vq) for several values of V Q (see fig. 4). This procedure was then 

not suitable, and (Vp^y was read into MAIN and remained unchanged for all of 
the values of Vq set by MAIN. 
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A list of the routines and a short description of their primary functions 
follows . 

MAIN determines what is to be computed (eigen-values or wave functions) and 
stores the values of the parameters to be used by the subroutines in computing 
the eigenvalues. 

Subroutine VOPFIX fixes the value of V M (x’) to be used in the computa- 
tions. 

Subroutine EVEIND sets limits on search for eigenvalues, stores the eigen- 
values found by other subroutines, and prints and/or punches the eigenvalues on 
sheets and/or cards. 

Subroutine ZERO locates changes in sign of the determinant as trial 

values for E are stepped, tests it for closeness to zero, and returns values 
of Dj^ and roots found to EVEIND. 

Subroutine CALC computes value of D^ for trial values of E given to it 
by ZERO and returns value of D-^ to ZERO. 


Details of Individual Routines 

MAIN provides needed flexibility in the computations. Desired values of 
N, w, h, a, and Vq as veil as the number of different values of Vq in a 
given run are fixed first. Then the following options are decided: 

VPSLCT determines whether the value of V^( x 1 ) will be computed by VORFIX 
or set in MAIN. 

EVSKIP determines whether eigenvalues or wave functions will be computed. 

EVODD determines whether even or odd solutions are to be used in the com- 
putations. If it is decided to compute eigenvalues, then DWRITE determines 
whether or not each computed value of will be printed out along with the 

trial value of E. 

Next, various quantities are read in that determine the starting trial 
value to be used for E, the amount by which subsequent trial values are to be 
stepped, and the quantities to be used to determine whether the computed value 
of Djj is close enough to zero to permit the corresponding trial value of E 
to be accepted as an eigenvalue. Control is then transferred to a DO-loop that 
sets values of v o after which the program is terminated. That portion of 
MAIN used to compute the wave functions will be described in appendix E. 

Subroutine VOPFIX, like MAIN, is used in computing wave functions as well 
as eigenvalues. Primarily, it sets specific values for Vq. Secondarily, it 
computes the corresponding values for H in kG as well as the number of the 
atoms in the chain at which weak breakdown and strong breakdown occur. If 
EVSKIP was not equal to 2 in MAIN, then VOPFIX will call EVEIND. 
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It should he noted that a dummy subroutine NORMAL must be included in the 
deck -when computing eigenvalues or the program will not run. 

Subroutine EVFIND sets the starting value of E at which the search for 

eigenvalues is to begin and the maximum value for E at which the search is to 

terminate. When an E and its corresponding DN value are received, both are 

stored, a new starting value is chosen to be given to ZERO, and the process is 

continued until all the eigenvalues in the desired range have been found. All 
of the eigenvalues and the corresponding values of DN are then printed. 

Those eigenvalues for which the DN f s are sufficiently close to zero are 
printed again separately. These latter eigenvalues are also punched out on 
IBM cards for use as input data in computing wave functions. 

Subroutine ZERO sends the starting value of E obtained from EVFIND to 
CALC, which sends back the corresponding value of DN. The value of E is 
then increased by an amount set in MAIN and called STEP. This new value of E 
is again sent to CALC, which again sends back the corresponding DN. The en- 
tire process is continued until the sign of the DN that is sent back to ZERO 
changes from the sign of the last DN sent. At this point, a linear interpo- 
lation procedure begins and continues until |DN| is smaller than an amount 
set in MAIN and called DTEST or until the amount by which the linear interpola- 
tion changes E is smaller than another quantity set in MAIN called PFfflCSN. 
When either of these two events occurs, the last values of E and DN are re- 
turned to EVFIND along with an indication (by means of sense lights) of which 
of the two events characterizes the particular pair of values of E and DN. 

Each trial value of E along with the corresponding value of DN will be 
printed or not in accordance with the value of DWRITE that was set in MAIN. 

Subroutine CALC utilizes the procedure indicated in appendix A to compute 
the value of DN when all of the parameters are fixed. Both even and odd so- 
lutions can be used in this subroutine. 


Input Data 

Units for the input data are the following: energy, electron volts; 

length, angstroms; magnetic field strength, kilogauss. 

The data are input from tape 5. The names of the quantities for which the 
data are used together with a description of the actual use are now listed. 

The subroutines in which these quantities are used follow each description. 

N number of atoms in positive half of chain (MAIN, VOPFIX, EVFIND, CALC) 

W width of well portion in one period of Vp(x ! ) (MAIN, CALC) 

H width of hill portion in one period of Vp(x') (MAIN, CALC) 

A distance of atomic separation (MAIN, VOPFIX, CALC) 

KI beginning value of cycle of values of Vq set in MAIN 
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KF 

KS 

VPSLCT 

VPCHS 

KIP 

KFP 

KSP 

EVSKIP 

EVODD 

DWRITE 

XSTART 

STEP 

PRECSN 

DTEST 

XTEST 

A 


largest value in cycle of values of Vq 

step in values of cycle of values of Vq 

1, has MAIN fix 2, has VOPFIX fix Vq (MAIN, VOPFIX) 

value of Vq given to MAIN (this card omitted if VPSLCT is 2) (MAIN, 
VOPFIX) 

beginning value of multiple of fractional value of (V' ) used in 

VOPFIX to set values of Vq (MAIN, VOPFIX) max 

final value of preceding description 

step in value of preceding description 

1, does not compute eigenvalues; 2, does compute eigenvalues (MAIN, 
VOPFIX) 

1, computes even eigenvalues; 2, computes odd eigenvalues (MAIN, CALC) 

1, values of DN printed in ZERO; 2, values of DN not printed in 
ZERO 

beginning trial value of E in a set of computations seeking zero 
value of DN (MAIN, EVFIND, ZERO) 

value by -which each successive trial value of E is stepped in pro- 
cess described previously (MAIN, EVFIND, ZERO) 

minimum interval of change from one trial value of E to another in- 
terpolated one that allows interpolation procedure to proceed (MAIN, 
ZERO) 

maximum absolute value of DN that permits the eigenvalue to be used 
in computing wave functions (MAIN, ZERO) 

minimum interval of change from one trial value of E to another in- 
terpolated one that allows ZERO to keep searching; if change is 
smaller than this and DN is less than DTEST, trial value is suit- 
able for use in computing wave function; if DN is still larger 
than DTEST, search stops but the last trial value of E is unsuit- 
able for use in wave function computations (MAIN, CALC) 

listing of the subroutines used in computing the eigenvalues follows. 
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SIBFTC WAIN DECK 

DIMENSION J DO ( 1 00 ) 

DOUBLE PRECISION XSTART, STEP, XMAX, PRECSN, DTEST, XTEST 
DOUBLE PRECISION U , H , A , VO , VPCHS 
INTEGER XDSIZE.EVSKIP.CWRITE, EVODD, VPSLCT 
COMMON OCR 

EQUIVALENCE I A, COM ( 1) J , (W,COM( 3) ) , (H,COM( 5) J , ( V0,C0M(7) ) , 

L (N«CCM(16) ) , ( XSTART , E S T AR T , COM ( ID) , ( VP SLC T , COM ( 330 ) ) , 

2 (DwRI TE.COM (796) ) , ( KIP, COM ( 13) ) , (KFP,COM( 14) ) , (KSP,COM{ 15 ) ) , 

3 ( N I , COR ( 17 ) ) , (NF,COM( 18) ) , (NS, COM ( 19) ) , ( EVODD , COM ( 32 9 ) ) , 

4(EvSKlP , COM (223)), (JEV »COM ( 224 ) ) , ( JDO, COM ( 225 ) } , ( KPLOT , COM { 574 ) ) , 
5(STEP,C0M( 797) ) , ( PR ECS N , COM ( 799 ) ) , ( XM AX , EMA X, COM { 00 1) ) , 

6(DTEST »C0M(803) ) , ( V PCHS , COM ( 33 1 ) ) , ( XDS I ZE , COM ( 00 7 ) ) , 
7(XTEST»C0M(811) ) , 1 LTM, COM( 8 14 ) ) , ( L TN, COM ( 8 l 5 ) ) , ( KN.COM ( 816) ) , 

8 ( KSX,CCM(817) ) , ( KS Y , COM ( 8 1 8 ) ) , ( FX , C OM ( 8 19 ) ) , ( DX , COM ( 020 ) ) , 

9 ( F Y , COM (821)), ( DY.COM ( 022) ) 

1 READ ( 5 , 7 1 ) N 

71 FORMAT ( 14 ) 

READ (5,72)W,H,A 

72 FORMAT ( 3D1C .3 ) 

READ 15.75JKI ,KF,KS 

75 FORMAT (314) 

C VPSLCT = 1,2 HAS MAIN.VOPFIX CHOOSE VOP 
READ (5,71) VPSLCT 
IF (VPSLCT. EQ. 2) GO TO 6 
READ (5,78) VPCHS 
78 FORMAT (D15.8) 

6 READ (5,75)KIP,KFP,KSP 

C EVSKIP = 1,2 MEANS WAVE FUNCTION, EV CALCULATED, RESPECTIVELY 
READ(5, 74) EVSKIP.JEV 
74 FORMAT (214) 

READ (5,71) EVCDD 

C E VC CD =1.2 GIVES EVEN, ODD SOL U T I ONS , R E S PEC T I VE L Y 

WRITE (6,62) (N, W,H, A, KI,KF,KS, KIP, KFP.KSP, EVODD, EVSKIP) 

62 FORMAT ( 1H K , 1 0 X , 2HN= I 2 , 2X , 2 HW= F4 . 2 , 2 X , 2 HH= F 4 . 2 , 2 X , 2HA = F 4 . 2 , 2 X , 
19HKI ,KF,KS=3I3, 14H, K I P , KF P , KS P= 3 I 3 , 8H , EVQDD=I2,9H, EVSKIP=I2) 

IF (EVSKIP. EQ.l) GO TC 4 
READ (5,71)CWRITE 
C DWRITE =1 WRITES DN 

READ (5,73)XSTART .STEP 

73 FORMAT (2D11. 4) 

READ (5,77) PRECSN, DTEST , XTEST 
77 FORMAT (3C11.4) 

WRITE (6, 63) (DWRITE, XSTART, STEP, PRECSN, DTEST, XTEST > 

63 FORMAT ( 1HK, 10X»7HDWRITE=I2»2X, 12HXSTART ,STEP=2D12.4, 
12X,19HPRECSN,DTEST,XTEST=3C12.4) 

GO TO 5 

4 READ ( 5 , 75 ) Nl » NF * NS 
READ ( 5 , 76 ) J DC 

76 FORMAT (70 1 1/3011) 

C KPLOT = 1,2 MEANS PLOT MADE, SKIPPED .RESPECTIVELY 
READ (5,71) KPLOT 
IF (KPLOT. EQ. 2 ) GO TO 5 
C NO OF POINTS IS = TO THE VALUE OF XDSIZE 
READ (5,71) XDSIZE 

C LTM SPECIFIES NUMBER LINE SPACES BETWEEN GRID LINES 
C L TN SPECIFIES NUMBER OF PRINT SPACES BETWEEN GRID LINES 
C KN IS THE NUMBER OF CURVES 

C EXPIKSX.KSY -6) TIMES F X , F Y OR TIMES DX.DY = ACTUAL STARTING VALUES 
C OR CHANGES IN GRID VALUES 

READ (5,51) LTW,LTN,KN,KSX,KSY 

51 FORMAT (514) 

C FX USED TO SPECIFY STARTING VALUE OF VERTICAL SCALE 

C DX USED TO SPECIFY CHANGE IN V/FRTICAL GRID VALUES EACH LINE SPACE 
C FY USED TO SPECIFY STARTING VALUE OF HORIZONTAL SCALE 
C DY USED TO SPECIFY CHANGE IN HORIZONTAL GRID VALUES EACH PRINT SPACE 
READ (5,52) F X , DX , F Y, DY 

52 FORMAT I4F8.3) 

5 IF (KI.EQ.O) GC TO 11 
DO 10 1= KI.KF.KS 

2 VO=I 

1C CALL VCPFIX 
GC TO 12 

11 VC = O-OCO 
CALL VCPFIX 

12 STOP 
END 



SUBROUTINE VOPFIX 

C GIVES OPTION OF NOT CALCULATING EV 

DOUBLE PRECISION Vi , H, A , VO , VOP , EV , VPCHS , XN , VP I NT , B VOP , C 
INTEGER VPSLCT, EVSKIP 
COMMON COM 

EQUIVALENCE { A,COM( 1) > , ( W , COM ( 3 ) ) » ( H, COM < 5 )), ( VO, COM ( 7)), 

1(N, COM (16 3 1 , ( VCP,COM( 325) ) , ( V PSLCT , COM ( 330 ) ) f ( VPCHS , COM { 33 1 ) ) , 
2(DWRITE,C0M(796) ) , ( K I P , COM ( 13 ) ) , ( KFP , COM { 14 ) ) , CKSP,C0M(15) ) , 

3 (EVSKIP, COM (223) ) , ( JEV , COM ( 224 ) ) 

IF (KIP.GT.O) GO TO 1 
VOP = 0 . DO 
WRITE (6,63) 

63 FORMAT ( l FK , 5X , 6HV0P=0 , 3X , 5HHM=0 , 3X , 2 1HN0 MAGNETIC BREAKDOWN/ LHK ) 
GO TO (16,17) .EVSKIP 

16 DO 18 J2 = 1 , J E V 
READ (5,71) EV 
CALL NCRMAL(EV) 

18 CONTINUE 
GO TO 15 

17 CALL EVFINC 
GO TO 15 

1 XN = N 

V P I NT - KFP 

BVCP = VC*CSQRT(U5D0)/(2.0D0*XN*VPINT) 

2 CH = 8.793984E-12 

DO 14 J= KIP, KF P , KS P 
C = J 

GO TO ( 19,3),VPSLCT 

19 VOP = VPCHS 
GO TO 20 

3 VCP = C*BVCP 

20 HM= SQRT( VCP/CHJ/A 

9 WRITE (6,62) VO, VOP, HM 

62 FORMAT ( 1HK ,5X»3HVO : =F5. 1»2X,4HV0P=F5.4» 16X, 3HHM-, 1PE9.3,/1HK) 

GO TO (11,13) .EVSKIP 

11 JEV -JEV 

DO 12 Jl= l , JE V 
READ ( 5 , 7 1 ) E V 
71 FORMAT ( D23 • 16 ) 

CALL NCRMAL(EV) 

12 CONTINUE 
GO TO 14 

13 CALL EVFINC 

14 CONTINUE 

15 RETURN 
END 



SUBROUTINE EVFIND 

DOUBLE PRECISION E , CN , ROOT , XST ART , XMAX , PREC SN , XBEG I N , ROOT I 
DOUBLE PRECISION E V , 6 V 1 , DE V , DEV 1 , W , H, A , VO, VOP , XTE S T 
D I PENSION EV( ICO) ,EVl ( 100) , DEV ( 100) , DEV U 100) 

COMMON COM 

EQUIVALENCE ( A.COMl 1) ) , (W,C0M(3) ) , (H,C0M(5) ) , < V0,CUM(7) ) , 

KNtCCM (16) ) , ( VCP.COMI 325) ) , ( ROOT • C0M{ 9 ) ) , ( XST ART , E START f COM 11 1 ) ) , 
2 ( NEV , CGM ( 2 0 ) ) , (EV.COM ( 21) ) , ( E , X , COM ( 32 7 ) ) , ( XTEST.COMt 811 ) ) , 

3 (STEP, COP ( 797) ) , ( PRECS N , COM ( 799 ) ) , ( XMAX , EMAX , COM ( 801 > ) , 

A (DTEST.COM (803) ) , ( ON , CCM ( 805 ) ) , ( XBEG IN , COM ( 809 ) ) 

DO 6 J=l f lCO 
EV(J) * 0. CO 
DEV(J) = 0.000 
DEVl(J) = 0 .000 

6 EV1IJ) = 0 • DO 

1 XN = N 

IF ( VO.EQ.O.ODO) GO TO 25 

2 BE MAX= XN*VO»SCRT( 1.5J/2.0 
IF (BEMAX-VO) 3,3,4 

3 EMAX = VC 4 0.5 
GO TO 5 

4 EMAX = BEMAX 4 0.5 
GO TO 5 

25 XMAX = 25.000 

5 XBEG I N * XSTART 
NEV = 0 

DC 15 J= 1,100 

7 CALL ZERO 

C SENSE LIGHT 1 CN IF NO ROOT FOUND FOR X .GE. EMAX 
CALL SLITET (1 .KOOOFX) 

GC TO ( 8 » 9 ) .KOOOFX 

8 NEV = NEV 4 1 
EVINEV) = C. DO 
DEV (NEV) = O.ODO 
EVl(J) = O.ODO 
PEVKJ) = O.ODO 

. 0 TO 16 

9 (J.GT.l) GO TO 11 
1C EV1(J) = RCOT 

DEV1IJ) = CN 

C SLITE 3 ON IF RCOT CANNOT PASS DN TEST 
CALL SLITET(3,K) 

GO TO (19,18) ,K 

C NEV = NUMBER OF EVS WHICH ARE CN TESTED 

18 NEV = NEV +1 
EV(NEV) = RCOT 
DEV(NEV) = DN 

19 R0CT1 = ROCT 
GO TO 14 

11 IF (ROCT • GE. EMAX ) GO TO 16 

12 I F ( DAB S (ROCT - ROOT l ) .GT. PRECSN) GO TO 10 
CALL SLITET(3,K) 

GO TO (13,22) , K 
22 IF (NEV.GT.O) GO TO 21 
NEV = NEV 4 1 
21 EV(NEV) * RCOT 
DEV(NEV) = CN 

13 XBEG I N = RCOT 4 STEP 
GO TO 7 

14 XBEGIN =* RCCT 

15 CONTINUE 

16 WRITE (6,64) VC.VOP 

64 FORMAT (1HK.10X, 3HV0=F4. 1 , 4X, 4HV0P=F8 . 4/ 1HK ) 

WRITE (6,61 ) (EV1( I) , DEVI II). 1=1, J) 

61 FORMAT (50X.9HALL ROOT S// 2 ( 28X , 2HE V , 28X , 2HDN ) / / ( 4030. 16 ) ) 

WRITE (6,62) (EV( I ) , DEV < I ) ,1-1, NEV) 

62 FORMAT ( 30X ,37HE IGENVALUES HAVING DN LESS THAN DTEST// 
12(28X,2HEV,28X,2HDN)// (4D30.16) ) 

DO 20 L=1,NEV 

20 WRITE (6,63) (EVIL) ,DEV(L) ,N,VO,VOP) 

63 FORMAT ( 1H$ , D23 . 16 , D 12 . 3, 2X , 2HN= I 2 , 2X, 3HV0=F4- l, 10X , 4HV0P=F 6. 3 ) 
RETURN 

END 



$ I6FTC ZERO LIST, REF, DECK 

SUBROUTINE ZERO 

DOUBLE PRECISION XSTART , STEP , XMAX, PRECSN, DTEST , XBEGI N, ROOT , DN 

DOUBLE PRECISION X, XI , X2 , T1 , T2 , DABS, E , XTEST 

DOUBLE PRECISION W,H,A,VO,VOP 

INTEGER DWRIIE 

COMMON COM 

EQUIVALENCE (A, COM! i) ) , ( W,COM( 3M , I H,COM( 5) ) , < VO , COM ( 7 ) ) , 

1 ( N , COM ( 16) ) , { VOP,CDM( 325) ) , ( ROOT, COM ( 9) ) , ( XSTART, ESTART , COM ( 1 1 ) > , 

2 ( STEP, COM { 79 7) ) , ( PREC SN , COM ( 799 ) ) , l XMAX, EMAX, COM( 801 ) ) , 

3 ( XTEST, COM ( 811) ) ,( DTEST , COM ( 803 ) ) , ( DN,COM ( 805 ) ), 

4( E, X, COM ( 327 ) ), ( DWR I TE , COM ( 796 ) ) , ( XBEGIN,COM( 809) ) 

X = X8EGIN 
CALL CALC 

GO TO ( 101,102) » DWR I TE 

101 WRITE (6,61) E, DN 

61 FORMAT! 10X,2(D28. 16) ) 

102 IF (DN) 1,2,3 

1 J = 1 
GO TO 4 

2 ROOT = X 

IFtDWRiTE.EQ.2) GO TO 15 

120 WRITE (6,63) ROOT , DN 

63 FORMAT! 10X,2(D28. 16) ,4HRU0T) 

121 GO TO 15 

3 J = 2 

4 XI = X 
T 1 = ON 

X = Xl+STEP 
CALL CALC 

I F ( DWR I TE • EQ* 2 ) GO TO 104 

103 WR ITE(6,61) E , DN 

104 IF (ON) 5,2,6 

5 GO TU (7,8) ,J 

6 GO TO ! 8 , 7 ) , J 

7 IF (X-XMAX) 4,14,14 

8 T2 = DN 
X2 =X 

9 X = ( DABS! T 1 ) »X2 + DABS ( T2 ) *X1 )/( DABS ( T1 i +DABS < T2 ) ) 

IF( X2-X.LE. PRECSN) GO TO 16 

10 CALL CALC 

I F ( DWR 1TE» EQ • 2 ) GO TO 106 
WR I TE ( 6, 61 ) E , DN 

106 IF (DN) 11,2,12 

11 GO TO ( 13, 8), J 

12 GO TO ! 8 , 1 3 ) , J 

13 XI = X 

T 1 = DN 

X = ( DABS! T1 ) *X2 + DAB S ( T2 ) *X 1 )/( DABS ( T 1 ) + DABS ( T2 ) > 

IF (X-Xl —PRECSN ) 17,17,10 

16 J2 = 2 

GO TO 18 

17 J 2- 1 

18 CALL CALC 

IF!DABS(DN).GT. DTEST) GO TO 20 

19 Jl=l 

GO TO 21 

20 J 1=2 

21 GO TO (2,22), J1 

22 GO TO (23,24i,J2 

23 IF ( X— X 1— X TEST )25,25, 10 

24 I F ( X2—X— XTEST ) 25 , 25 , 10 

25 ROOT = X 

GO TO ( 107,108) , DWR I TE 

107 WRITE (6,62) E, DN ** 

62 FORMAT ( 10X , 2 ( D28 . 16 ) , 1 2HDN TOO LARGE) U 

108 CALL SLITE (3) 

GO TO 15 

14 CALL SLITE ( 1 ) 

15 RETURN 
END 


56 



SUBROUTINE L AC C 
INTEGER EVCLC 
SUBROUTINE LALC 
INTEGER EVCLC 

DOUBLE PRECISION «, H, A , VU, VOP 

DOUBLE PRECISION UN , E , DA B S , DSQR T , GO w , WGO W , RGO W , R BGOW f X I f X J , XN 
DOUBLE PREClSlCNBGJH,BGIrf,AUGJH,ABGU,GJH,GIW,RGl*,RGJH,HBGJH 
DOUBLE PRECISION LAPP , D£XP , THGJM * EXPN , C SHF 
DOUBLE PRECISION SNHF , CSF , DCOS , SNF , DSIN 
DOUBLE PRECISION R BG l H , XN, BGNH, A BGNH , GNH, E T A f C S , T WG I W, T WGNH 
CUMMCN COM 

EQUIVALENCE (-A,COM( l ) i , ( H» COM I 3) ) , ( H , COM ( 5 ) I , ( VO, COM ( 7)1, 
l IN, COM | 16)), (* VOP , COM ( 325 M , ( ROOT , COM! 9 I ) , I X S T AR T , E S T AR T , COM { III), 
2(*E*X,C(JM( 32 7 M * I EVODD.COMI 329) ), I ON .COM I 805) ) 

EQU I VALENCE ( RGOw,RGI w) , ( RBGOW, RBG I W ) 

C SUBROUTINE SYMUUL6 RG I W, RBG I H, RGJH, RBGJH CORRESPOND TO TEXT 
C SYMBOLS SN h » SN*-B AR , SNH , SNH- BAR .RESPECTIVELY 

701 ETA =0-512335100 

702 GO h = ETAtObCRTIE) 

703 WG0W***G0* 

GO TO ( 70<t, 705) , EVOOD 
70* RGO W *- DCOSI WGO W ) 

RBGOm *GOm«USIN( WGOM) 

GO TC 200 

705 RBGOW « -OCGSlHGO*) 

IF (V.GOW -Em- 0.00) GO TO 706 
RGO W -= -DSIN(’WG0W)/G0W 
GO TO 2C0 

706 RGO W * -rf 

200 DO 300 1=1, N 
J« 1-1 

X I * I 
XJ =J 

BGJH= E- X J # »2 * VOP- VO 
BG l R * E- X i • * 2 • VOP 
A0GIw=OABSIdGIW) 

ABGJH=CABSIBGJH) 

GJH = ETAfCbCRT ( ABGJH) 

GIW * ETA*OSGRT ( AdG 1 w ) 

TMGIW * 2.00 » h * U I W 
THGJH « 2. CO »H« GJH 
IF ( BG JH ) 1,2,3 

1 E X P P= DEXPI IHGJH) 

EXPN = 0 E XP I — T HG JH ) 

CSHF = (EXPP* E X PN ) / 2 • 00 
SNhF = ( EXPP- E XPN ) / 2 . 00 

RGJH * - I KG I W »CSHF ♦ R BG I W* SNHF / G JH ) 

RBGJH = -I HGtW«GJH*SNHF *RBGIW*CSHF) 

GO TC 210 

2 RGJH =- I RGIW+2.U0«H«RBGI H ) 

RBGJH *- RBG I W 

GO TO 2 L 0 

3 CSF = CC0S4 IHGJH) 

SNF = DSINI Tt-GJH)- 

RG JH * -I RG I W *C SF * RBGIW»SNF/GJH) 

RBGJH = RGlRtGJH » SNF - RBGIW*CSF 
210 IF(BCIW) 201,202,203 

201 EXPP = CEXPI IWGIW) 

EXPN * UEXPI-TWGI*) 

CSHF = (EXPP* E XPN 1/2.00 
SNHF - (EXPP- E XPN ) / 2 . DO 

RG I K = -IRGJH*CSHF ♦ RBGJH*SNHF/GIH) 

RBG I W = — ( RGJH»GIW« SNHF *RBGJ H» C SHF ) 

GO TC 300 

202 RG I W =— I RGJH*2.D0«W*RBGJH) 

RBG I H *- RBO J H 

GO TC 300 

203 CSF = DCUSI IWGI W ) 

SNF = CSINIThGlH) 

RG I W • - 1 R GJH* C SF ♦ R BG J H* S NF / G I W ) 

RBGIH = RGJH *GI**SNF -KBGJH*CSF 

300 COM I NU t 
XN -N 

BGNH «= £— XN*»2* VOP— VO 
ABGNH= DABS ( BGNH ) 

GNH = ETA*DbCRT( ABGNH) 

THGNH * 2.00 •H»GNH 
I F I BGNH ) 301,302,303 

301 E X PP = UEXP1 IHGNH) 

EXPN «= DEXPl-THGNh) 

CSHF = (EXPP* E XPN ) / 2 . DO 
SNHF ■= (EXPP- E XPN ) / 2 . DO 

311 ON = — (RGlRfCSHF* R BG I W* SNHF / GNH ) 

GO TC 6C0 

302 DN«-I RG IW*2.C0*H«KBG I W) 

GO TL 600 

303 CSF = CCDs I IHGNH) 

SNF = CSINl THGNH) 

313 UN = -(RGIwtCSF • RBG I W*SNF/GNH) 

600 RETURN 
F ND 



APPENDIX E 


FORTRAN IV PROGRAM FOR COMPUTATIONS OF WAVE FUNCTIONS 
General Description 

Like the eigenvalue program, this program operates vith the Lewis 7094 
monitor system. It consists of the same MAIN listed in appendix D plus eight 
subroutines. The following is a description of the primary functions of the 
routines : 

MAIN determines what is to be computed (eigenvalues or wave functions) and 
stores the values of the parameters to be used by the subroutines in computing 
the wave functions. 

Subroutine VOPFIX fixes the value of V-^(x l ) to be used in the computa- 
tions and feeds eigenvalues to subroutine NORMAL. 

Subroutine NORMAL computes and stores ratios of coefficients of sin- and 
cos-like terms in each interval (the A n , s and B n T s in eqs. (45) to (47)) to 
the arbitrary coefficient in the 0 th well. (These quantities are the a n 's 
and (3 n T s of eqs. (C15) to (C18).) It also normalizes the wave function by 
finding the value of the (or Ag) that will make 

r°° 

/ |A(x')| 2 dx' = 1 

•/-00 

in accordance with the procedure in appendix C. 

Subroutine WFC computes the normalized A^'s and B n T s and sends them to 
other subroutines to be used in computing actual values of the wave function in 
the entire range x T = 0 to (N + l)a - w. The other subroutines return the 
wave function values to WFC, which prints them all out and determines whether 
or not a plot should also be made. 

Subroutine WFO computes values of A^(x f ), A^(x ! ), |A^(x t )| 2 j and 
|Aq(x t )| and returns these values to WFC. 

2 

Subroutine WFXIW computes values of A^(x f ) and |A^(x ! )| for n > 0 
and returns these values to WFC. 

Subroutine WFXIH computes values of |aJJ(x')| and |aJ(x')| 2 for n>0 
(including n = N) and returns these values to WFC. 

Subroutine PLOTPX sets up PLOTMY. 

Subroutine PLOTMY furnishes a plot of A(x') and |A(x')| 2 against x'. 
(This subroutine is part of the library tape of the Lewis Monitor System.) 
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Details of Individual Routines 


In the following descriptions only those parts of MAIN and VOPFIX pertain- 
ing to computations of the wave functions are included. 

In MAIN the parameters N, v ; h, a, and Vq as well as the choice of how 
V^(x I ) is to he made are fed in as described in appendix D. EVSK1P must be set 
equal to 1 so that eigenvalues will not be computed. In addition, JEV must be 
fixed, which gives the total number of eigenvalue data cards to be used in the 
run. 


Next, the number of subintervals in each region of constant potential is 
fixed. Provision is made for computing wave functions only in certain regions 
if desired. Then a’ decision is made as to whether a plot is to be obtained in 
addition to the printed values for A(x'). If a plot is asked for, various 
parameters needed by PLOTMY are then read in. Control is then transferred to 
the same DO-loop that sets Vq as described in appendix D at the end of which 
the program is terminated. 

Subroutine VOPFIX will determine the value of V M (x T ) to be used in subse- 
quent subroutines as described in appendix D. However, if EVSKIP was set equal 
to 1 in MAIN, then instead of calling EVFIND, control will be transferred to an 
interval DO-loop that will read eigenvalues stored in MAIN one at a time and 
call NORMAL. When this procedure has been followed JEV times (see MAIN), the 
value of Vjy^( x * ) in the external DO-loop is stepped up and the inner DO-loop 
cycle repeated. When the outer DO-loop is completed, control is returned to 
MAIN. 


It should be noted that a dummy subroutine EVEIND must be included in the 
deck when computing wave functions or the program will not run. 

Subroutine NORMAL computes A^/B^, B^/B^, A^/b^, and b£/B™ (see eqs. 

(C13) to (C18)) from the matching equations (A4) , (A5), (A7), (A8) , (A9), 

(AiO), and (All). As explained in the text, if DN is small enough for the 
eigenvalue being used, the results will be the same as if equation (A12) had 
been used instead of equation (All). The routine is set up in such a manner 
that individual quantities needed for the normalization as given by equations 
(C21), (C22), (C23), (C24), and (C9) are computed and stored along with each 
ratio. After all of the ratios have been computed and stored and the normaliz- 
ation of the arbitrary coefficient has been completed, WPC is called. 

Subroutine WPC normalizes and stores the unnormalized coefficients com- 
puted in NORMAL. Then A^, B^, and are computed and WPO is called. Next, 

a DO-loop is entered in which the following procedure is carried out: 

First, a test is made to determine whether A(x f ) was desired for the re- 
gion being considered. If the data in MAIN indicated that no X(x*) was asked 
for the region in question, then the same test is made for the next region one 
atomic distance farther along the positive* half of the chain. If the test in- 
dicates that X(x’) was asked for in the region, then A^, B^, and are 

set up and WFXIW is called to compute ^(x 1 ) in the well portion of the region. 
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When WFXIW returns control to WK) ; A|*, B^, and are set up and WFXIH is 

called to compute A(x') in the hill protion of the region. Following this, 
the test is applied to the next region. 

After all of the regions in the chain have been tested and values of 
A(x’) computed for all of the region for which these values are requested, the 
values are printed out. 

Then if a plot is desired (KPLOT = l) , PLOTFX is called. If no plot is 
desired, control is returned to NORMAL. 

The special subroutine WFO was written to compute Aq(x t ) and Aq(x* ) be- 
cause the beginning and end of the cycle have somewhat special behavior. 

These quantities and their absolute squares are computed by using equa- 
tions (24) and (25). Control is then returned to WFC. 

Subroutine WFXIW computes the quantities A^(x ! ) and [A^x 1 ^ (for 
0 < n <; N) using equation (38). Control is returned to WFC. 

Subroutine WFXIH computes the quantities A^(x T ) and jA^(x T )| 2 tor 
0 < n < N by using equation (39). When a test indicates that n = N, then 
A^(x ! ) and |A^(x T )[ are computed by using equation (40). After each set of 
computations control is returned to WFC. 

Subroutine PLOTFX(EV) arranges the values of x T , A n (x ! ), and |A n (x T )|^ 

in arrays so that they can be plotted properly by PLOTMY. All of the input 
control data for this subroutine is read in by MAIN. The eigenvalue is an 
argument of this subroutine and is read in by VOPFIX. Actually, the eigenvalue 
is only required for the legend of the plot. 

Subroutine PLOTMY is part of the library tape of the Lewis monitoring sys- 
tem. It is the routine that actually plots A(x T ) and |A(x f )|^ against x T . 
The main feature of this routine that must be taken into account is that it 
places the ordinate across the top of the page with zero to the left and the 
abscissa down the page with the lowest value at the top. While this feature is 
somewhat inconvenient in many cases, it has the advantage of permitting the 
case of an abscissa of arbitrary length since several pages can be covered con- 
tinuously. Its use is fully described in reference 16. 


Input Data 

Units for the input data are the following: energy, electron volts; 

length , angstroms . 

The data are input from tape 5. The names of the quantities for which the 
data are used along with a description of the quantity are now listed. The 
subroutines in which these quantities are used follow each description. 

N number of atoms in positive half of chain (MAIN, VOPFIX, NORMAL, WFC, 

WFXIH, PLOTFX) 
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w 


width of well portion in one period of Vp(x') (MAIN, NORMAL, WO, 
WXIW) 


H width of hill portion in one period of Vp(x j ) (MAIN, NORMAL, WO, 

WXXH) 

A distance of atomic separation (MAIN, VOPFIX, NORMAL, WXIW, WXIH) 

KI see appendix D 

KF see appendix D 

KS see appendix D 

VPSLCT see appendix D 

VPCHS see appendix D 

KIP see appendix L 

KFP see appendix D 

KSP see appendix D 

IVSKIP see appendix L 

EVODD see appendix D 

NI beginning multiple of fraction of subdivision of individual region 

used for computing wave functions 

NF l/NF, fraction of subdivision of individual region; NF, final multiple 

of this subdivision 

NS stepping interval of multiple of subdivision 

JDO dimensional quantity by means of which the decision to compute Tk(x') 

for a region is made; each region in sequence is characterized by 
the location of the number on the data card - the first number on 
the card corresponds to the 0^ region; if the column on the 

data card is blank, no computation of Aj(x') will be made; if the 
column contains a 1, the computation is made (MAIN, WC) 

KPLOT 1, plot of A(x r ) will be made; 2, no plot of A(x’) will be made 

XDSIZE number of points on one curve in plot 

LTM number of line spaces between grid lines on plot 

LTN number of print spaces between grid lines on plot 

KN number of curves on plot 
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KSX scaling parameter for x-scale (runs up and down the page); FX and DX -will 
he multiplied hy loKSX-6 

FX quanti J ' r used to specify starting value on vertical scale; actual start- 
ing value , FX times lo^SX-6 

DX quantity used to specify change in value in vertical scale one line 
space; actual change, DX times io^SX-6 

KSY scaling parameter for y-scale (runs across page); FY and DY ■will he mul- 
tiplied hy loKSY-6 

FY quantity used to specify starting value on horizontal scale; actual 
starting value, FY times l()KSY-6 

DY quantity used to specify change in horizontal scale in one print space; 
actual change, DY times e^SY-6 

A listing of the subroutines used in computing the wave function follows 

(MAIN and VOFFIX are listed in appendix D) . 
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SUBROUTINE NURMAL(EV) 

C FULL 00U3L £ PRECSN. UNLY STORES FOR 50 ATOMS NOW. 

DIMENSION AKWI50I , BKW( 50) , AKH { 50) , BKH< 50) , 
lGKrft 50) ,GKH( 50) , J W ( 50 ) f J H ( 50 ) ,RMLZ( 105) 

DOUBLE PRECISION W, H, VO , VOP , BOW, G I W, GJH, AKH , BKW, AKH, BKH, GK W , GKH 
DOUBLE PRECISION OAbS » OSQR T * E V , EX , BG JH , ABG J H , BG I W , ABG I W , X I , R ML 


DOUBLE 

PRECISION 

SNF 

f 

OSIN 

» 

TWGI W 

• 

EXPP 


DEXP 

DOUBLE 

PRECISION 

FHGJH 


EXPN 

t 

SNHF 

» 

CNHF 


CNF 

DOUBLE 

PRECISION 

OCQS 

» 

WGIW 

* 

EXPPJH 

t 

THGJH 


EXPNJH 

DOUBLE 

PRECISION 

SNHF J H 

t 

CNHF JH 

f 

EXPPIW 

» 

FWG I W 


EXPNI w 

DOUBLE 

PRECISION 

SNHF 1 W 

• 

CNHF I W 

» 

SNF JH 

♦ 

CNF JH 


SNF I W 

DOUBLE 

PRECISION 

CNF 1 w. 

RMLZ, ETA 







DOUBLE 

PRECISION 

AANMOH, 

BBNMOH, ABNMOH, A ANM I W 

, BBNM I W , 

ABNMIW 


DOUBLE PRECISION A ANM J H , BBNM J H , ABNMJ H 
INI EGER EVOOD 
COMMON COM 

EQUIVALENCE ( A* COM( 1) ) , ( W,COM( 3) ) , (H,COM( 5) ) , ( VO,COM( 7) ) , 

I ( VOP, CUM ( 32 5) ) , I BOW, COM I 950 3) ),( AKW, COM! 9505 ) ) , ( BK W , C DM ( 960 5 ) ) , 

21 AKH, CUM (9705) ) , ( BKH , COM ( 9 805 > ) , ( GKW, CUM ( 9905 ) ) , ( GKH , COM ( 1 00 0 5 ) ) , 
3 ( N , CUM ( 16) ) , ( JH.COM ( 10105) ) , (JW.COM ( 10205) ), ( E VODD, COM ( 329 ) ) 

1 tX = EV 

ETA = 0.512335100 
AK W { 1) = O.OUO 

BKW< 1 ) = O.OUO 
GIN = ETA*USQRT ( EX) 

GKW( 1) = GIW 
WGIW = W*G I W 
rWGIW = 2 • DO* WG I W 
BG J H = EX-VO 
A BG JH = DABS ( BG J H ) 

GJH = tTA*DSuRT ( ABGJH) 

GKH ( 1) = GJH 

T HG JH = 2 . DO *H* G J H 
FHGJH = A . DO *H* G J H 
Jrf( 1) = 3 

CCUMPU I E COEFFICIENT OF BOW OR AOW IN RMLZ 
SNF = DSlN(TwGlW) 

C TERM IN NORMALIZATION INVOLVING BOW OR AOW 
GO TO l 101, 102) , EVOUD 

101 RMLZ(l) = W*( l.DO+SNF/TWGIW) /2.D0 
C COMPUTE BOH/BOW 

BKH ( 1 ) = DCOS ( WG I W ) 

GO TO 103 

102 RMLZ (1) = W* ( 1 * DO -SNF/TWGIW)/ { 2.U0*GIW**2) 

C COMPUTE BOH/AOW 

BKH ( 1 ) = DSINl WGIW) /GIW 

C COMPUTE CUEFFS OF AOH* * 2 , BOH# * 2 , AND A OH* BOH IN RMLZ 

103 IF (bGJH) 2,3,4 

2 JH( 1) = 1 

EXPP - DEXP { FHGJH ) 

EXPN = DEXP ( -FHGJH ) 

SNHF = ( EXPP-EXPN)/2.D0 
CNHF = ( EXPP+EXPN ) /2. DO 
C ( 1 /H ) * COEFF OF (A0H)**2 IN RMLZ 
AANMOH a SNHF /FHGJH- 1 . DO 
C (1/H)* COEFF OF oOH* * 2 IN RMLZ 
BBNMOH = SNHF/FHGJH+l.DO 
C (1/2H)* COEFF UF AOH* BOH IN RMLZ 
ABNMOH = ( CNHF- 1 .DO)/FHGJH 
GO TO 5 


3 Jri( 1) = 2 

AANIMOH = 2.DO*( 2,D0*H) **2/3. DO 
ttdNMOH = 2. DU 
ABNMOH = 2 - DO*H 
GO TO & 

4 JH( 1) = 3 

SNF = L)S IN { FHGJ H ) 

CNF = UCOS(FHGJH) 

AANMOH = 1 • DO— SNF/ FHGJ H 
B8NM0H - I . DO+SNF/FHGJH 
ABNMOH = ( 1 .OO-CNF ) /FHGJH 

5 GO TO ( 106,107) ,EVOL>D 
C CGMPUTt AOH/BOW 

106 SNF = USIN(WGIW) 

AKH ( 1 ) = -G I W*SNF 
GO TO 108 

C COMPUTE AOH/AOW 

107 AKH ( 1 ) = DCOMWGIW) 

108 JHK = J H ( 1 > 

GO TO (6, 7, 6), JHK 

6 AKH ( 1 ) = AKH( 1 ) /GJH 

C TERM IN RMLZ OUE TO ZERUTH HILL 

7 RMLZ l 2 ) = H* ( AANMOH*AKH( 1)**2+BBNM0H*BKH( 1)**2 
1 + 2 • DO * ABNMOH* AKH { 1 ) *BKH ( 1 ) ) 

DO 40 I = 1 , N 
XI = I 
II = 1+1 

BGId = EX-X I ** 2-* VOP 
ABGIrf = DABS(BGIW) 

GI «i = tTA*OSgRT( ABGIW) 

GKW(Il) =G I W 
TwGIW = 2 • 00 * vV* G I W 
FrtGlW = 4*00*W*GIW 
GO TQ< B,9, 10) , JHK 
C COMPUTE Blw/bOW 

8 EXPPJH = DEXP(THGJH) 
tXPNJH = OEXP(-THGJH) 

^NHFJH = ( EXPPJH-EXPNJH) /2.00 
LNHFJH = ( EXPPJH+EXPNJH)/2.D0 
BKW(Il) = AKH( I ) *SNHFJH+8KH( I )*CNHFJH 

AK W ( II) = GKH( I )*{ AKH(I)*CNHFJH+BKH( I )*SNHFJH) 
GO (0 11 

9 liKWtll) = 2.D0*H*AKH( I)+BKHl I) 

AKW (ID = AKH( I) 

GO TO 11 

10 SNFJH = OSIN(THGJH) 

CNFJH = DCUSITHGJH) 

BKW(tl) = AKrll I ) * SNF J H+6KH ( I ) *CNF J H 
AKW(Il) = GKH( I )*( AKH( I ) *CNFJH-BKH( I )*SNFJH) 

11 IF (BGirt) 12,13,14 

12 J W ( II) = 1 

C COMPUTE COtFFS OF A I rt* * 2 , B I w* * 2 , A I w* B I vi IN RMLZ 
LXPPIW = OEXP(FWGIW) 

EXPNIW = DEXP(-F^GIW) 

^NHFIW = { EXPPI W-EXPNI W) /2.D0 
LNHFIrt = 1 EXPPI W+EXPNIW)/2,00 
C ( 1 / w ) *COEFF OF A I W* * 2 

AANMlri = SNHF I W/ FWG 1 rt- 1 . DO 
C ( 1/W) *CU t FF OF 131 W**2 

bb'JBIw = SNHFI'J/FWGIW+i.DO 
C (1/2W)*C0EFF OF AlW*i3lW 

AdNMIW = ( CNHF I W— 1 • DO ) / F WG I W 
GO TO 15 
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13 J w ( II) = 2 

AAJMlW = 2.U0* ( 2*00*W ) **2/ 3 . DO 
ObNMlW = 2. DU 
A JNMIW = 2 • D 0 * W 
GU TO 15 

14 J w ( 11) = 3 

S.MFIW = USIN(FWGIW) 

CnFIW = UCO 5 ( F WG I W ) 

AAmMIW = 1 • OO-SNF I W/ F WG I W 
DBNHIW = I.DO+SNFIW/FWGIW 
ABNMIW = ( 1 .uO-CNFI W)/FWGIW 

15 J W K = J w ( II) 

C COMPUTE A I w/ BOW 

GO TO ( 20, 2JL, 20 ) , JWK 

20 AKW ( ID = AK W ( I 1 ) / G I Inl 

G TERM IN RMLZ DUE TO I TH WELL 

21 12 = 2 * 1+1 

KML Z ( 12) = W*( AANMIW*AKW( 1 1 ) **2+BBNMI W*BKW{ III **2 
1 + 2 • 00 * A BlNM 1 W*AK W( 1 1 ) *BKW{ I 1 ) ) 

C COMPUTt QUANTITIES FOR ITH HILL 

22 BGJH - BGIW-VO 

A BGJH = DABS(BGJH) 

G J rl = ETA*D$qRT ( ABGJH) 

GKH ( ID = GJH 
THGJH = 2 • DO*H* G J H 
FHbJH = 4 • DO * H* G J H 

C COMPUTE CUtFF OF AJ H* *2 , B JH** 2 , AND A J h*B J H IN RMLZ 
IF (BGJH) 23,24,25 

23 JH( ID = 1 

C EXPPJri AND EXPNJH NOT SAME AS THOSE USED BETWttN 8 AND 9 
EXPPJH = D6XP ( FHGJH) 

EXPNJH = DEXP(-FHGJH) 

SNHFJH = ( EXPPJH-EXPNJH) /2.D0 
CNHFJH = ( EXPPJH+EXPNJH) /2. 00 
L { 1/H)*CUEFF OF A J H **2 

28 AANMJH = SNHF JH/FHGJH- 1 . DO 
C ( 1/H ) *CQEFF OF BJH **2 

BBNMJH = SNHF JH/FHGJH+ 1 . DO 
ABNMJH = ( CNHFJH-1. DO) /FHGJH 
GU TO 33 

24 JH{ I 1) =2 

30 AANMJH = 2. DO* ( 2 . QO*H ) **2/3 . DO 
BBNMJH = 2. DO 
ABNMJH = 2 • DO*H 
bU TO 33 

25 J H ( ID =3 

SNFJH = DS I N ( FHGJH ) 

CNF JH = DCUi(FHGJH) 

32 AANMJH = l.DO-SNFJH/FHGJH 
BBNMJH = I .DO+SNF JH/FHGJH 
ABNMJH = .( 1 . dO-CNFJH) /FHGJH 

33 JHK = JH ( I l ) 

bO TO (27,29,31) , JWK 


65 


J 



C EXPPIW AMO EXPNIW DIFFERENT FROM THDSE BETWEEN 12 AMD 13 
27 EXPPIW = DEXP(TWGIW) 

EXPNIW = DFXP(-TWGIW) 

SMHFIW = (EXPPIW— EXPNIWI/2 . DO 
CNHFIW = (EXPPIW+EXPNIW)/2.D0 
IF ( I-N.GE.O) GO TO 41 
C COMPUTE BJH/BOW 

BKHUD = AK W ( II) *SNHF I W + BKW ( I i ) *CNHF 1 W 
C COMPUTE AJH/BOW 

AKH(Ii) = GIW*( AKW( ID*CNHFIW+BKW( I1MSNHF1W) 

GO TO 26 

liKH(Il) = 2.00*W*AKW( tl)+BKW{ I 1) 

29 IF ( I-N.GE.O) GO TO 41 
AKH ( ID = AK W ( I 1 ) 

GO TO 26 

31 SNFIW = DSIN(TWGIW) 

CNFIW = DCOSfTWGIW) 

IF ( I-N.GE.O) GO TO 41 

dKH ( II) = AKWU1)*SNFIW+BKW{ U)*CNFIW 

AKHCIl) = GIWM AKW( 1 1 ) *CNF I W-BKW ( I1)*SNFIW) 

26 GO TO ( 38,39,38) » JHK 

38 AKH ( ID = AKH ( ID/GJH 

C TERM IN RMIZ DUE TO ITH HILL 

39 13 - 12+1 

RMLZII3) = H*(AANMJH*AKH( I 1 ) * *2+BBNMJH*BKH { ID **2 
l+2.D0*A6NMJH*AKH( 1 1 ) *BKH( ID) 

40 CONTINUE 

C COMPUTE QUANT I TES FOR LAST HILL 

41 N 1 - N+l 

3KH( N1 ) = 0.00 

GO TO 142,43,44) ,JWK 

42 AKH ( N l ) = AKwtNl )*SNHFIW+8KW(NL)*CNHFIW 
GO TO 45 

43 AKH ( N 1) = 2.00*W*AKW( N1 ) +8KW( N1 ) 

GO TO 45 

44 AKH ( N l ) « AKrt{ND*$NFIW + BKW{ND*CNFlW 

45 IF ( 6G JH ) 46,47,48 

46 EXPPJH = OEXP ( T HGJH) 

EXPNJH - OEXP ( — THGJH) 

SNHFJH = ( EXPPJH-EXPNJH)/2.00 
AKH(Nl) = -AKH{ ND/6NHFJH 

C { 1 / H ) *C0£ FF OF AKH(N+1)**2 IN RMLZ 
EXPPJH = OEXP { FHG JH ) 
tXPNJH = OEXP ( -FHGJH) 

SNHFJH = ( EXPPJH-EXPNJH)/2.D0 
AANMJH = SNHFJH/FHGJH-l.DO 
GO TO 49 

47 AKH(Nl) = -AKH(N1)/(2.00*H) 

AANMJH = 2.00*( 2.D0+H) +*2/3. DO 
GO TO 49 

48 SNFJH = L)SIN(THGJH) 

AKH(Nl) = -AKH( N 1) / SNFJH 
SNFJH = L)S I N ( FHGJH ) 

AANMJH = l.DO-SNFJH/FHGJH 

C LAST ENTRY IN RMLZ 

49 N2 =2*N1 

RMLZIN2) = H*AANMJH*AKHIN1)**2 
RML = 0.00 
DO 50 J = 1 , N 2 

50 RML = RML + RMLZ(J) 

51 dOW = OSQRT { 1.00/ ( 2.D0*RML) ) 

CALL WFC { E V ) 

RETURN 

END 
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SUBROUTINE WFCtEV) 

DOUBLE PRECISION AK W , BKW , AKH, BKH , GK W , GKH, E V , VO, VOP 
DIMENSION JDO< LOO) , AKW (50) , BKW ( 50) , AKH( 50) ,BKH< 50) , 

IGKto ( 50 ) , GKH ( 50 ) , XP ( 220 ) , WF { 220 ) , WFSQ< 220) 

INTEGER EVCDD 
COMMON COM 

EQUIVALENCE i A,COM( 1) ) , (W. COMO) ), (H,C0M(5) ),( VQ,COM{ 7) ) t 
KN.COM (16) ) , (VC P, COM (325) ) , ( ROOT , COM ( 9 ) ),{ XSTAR T , E START f COM ( 1 1 ) ) , 
TCNI.COM (17) ). (NF.COMCIBK, (NS, COM { 19) ) . C JDQ .COM ( 225 > ) , 

2 t KPLOT .COM (574) ) . UK1.C0MC575) ). ( XP , COM ( 576 ) ) , (WF.COM (896) ) , 

3( WFSQ.COMC 1116) ) . (BOW, COM ( 9503) ) , ( AKW , COM ( 9505 ) ) , ( BKW , COM ( 9605 ) ) , 

4 (AKH,CCM(9705 )) , (BKH,CCM( 9805 ) ) , (GKW,COM( 9905) ) , ( GKH, COM ( 10005) ) , 

5 (EV0DD»CCM(329) ) 

GO TO (l.iCD.EVODO 

1 BKWU) = BCW 
GO TO 102 

10 1 AKW(i) = BCW 

102 N 1 = N+l 
DO 2 J= 1 , N 
Jl= J+l 

AKW(Jl) = BOW* AKW ( J 1 ) 

BKW(Jl) = B0W*BKW ( J 1 ) 

AKH(J) = BO W* AKH ( J ) 

2 BKH(J) = BOW* BKH(J) 

AKH(Nl) = BOW* AKH ( N 1 ) 

BKH(Nl) = O.ODO 

C ALL CCEFFS NOW NORMALIZED AND STORED 
IK1 = 0 

IF ( JDC(l) ) 4,4,3 

3 IK1=IK1 +1 
ANH = AKH (1) 

BNH = BKH I l ) 

GOW = GKW ( 1 ) 

GOH = GKH ( 1 ) 

CALL WFO( ANH, BNH, GOW, GOH) 

4 JO = NF/NS 
DO 10 1-2 * M 

IF ( JDC ( I ) ) 10,10,5 

5 IF ( I— N 1 ) 7,6,6 

6 CALL SLITE (2) 

7 IK1 = I K 1 +1 
ANW = AKW ( I ) 

BN W = BKW ( I ) 

G I W = GKW ( I ) 

CALL WFXlWtANW, BNW , G iVi , 1 ) 

IF U-M) 8,11,11 

8 ANH = A KH ( I ) 

BNH = BKH ( I ) 

GJH = GKH { I ) 

CALL WFXIH ( ANH, BNH, GJH, I ) 

1C CONTINUE 

11 CALL SLITET (2.KOOOFX) 

GC TO ( 12, 13) * KOOOF X 

12 ANH = AKH ( M) 

BNH = 0.0 
GJH = GKH(M) 

CALL SLITE (2) 

CALL WFXIH(ANH,BNH,GJH, I ) 

13 WRITE (6,61) VC ,V0P,EV,B0W, EVODD 

61 FORMAT ( 1HK»10X,3HVO=F5. 1 , 5X , 4HV0P=F 1 5 . 9, 5X , 3HE V= 1PD2 3 . 1 5 , 5X , 
14HB0W=1PE 15.7,5X,6HEV0CD= 12/// 

22 ( 4X, 3HXPW,6X ,4HWF I W, 4X,6HWFI WSQ , 4X , 3HXPH, 7X, 4HWF I H , 4X , 6HWF I HSQ ) ) 
NPLCT = N/2 +1 
DO 14 I A= l.NPLOT 
K= 40 * ( I A— 1 ) +1 
L-K + 9 

WRITE (6,63) 

63 FORMAT (1HK) 

WRITE (6,62) ( (XP( J) , WF ( J ) , WFSQ ( J ) , XP ( J ♦ 10 ) , WF ( J+l 0 ) , WFSQ{ J+10) , 

1XP( J+20) , WF ( J + 20 ) ,WFSQ( J+20),XPl J + 30) ,WF( J+30), WFSQ( J+30) } 

2, J=K, L ) 

62 FORMAT ( 4 ( C PF8 . 3 , 1 PE 10 .2, 1 PE 1 0. 2 ) ) 

14 CONTINUE 

GO TO (15, 16), KPLOT 

15 CALL PLOTF X ( E V ) 

16 RETURN 
END 



SUBROUTINE MFOC AHN, BHN , GOW , GOH) 

DIMENSION xe('2201,WFC220), WFSQC 220) ,JH( 1001 
INTEGER EVGGD 
COMMON COM 

EQUIVALENCE eM*C0M{3) ) , CH,C0H(5) ) , INI, COMC 171) , CNF, COMC 18)1, 
lC'NS, COM { 1911 *I XR,C0H<576)) , CWF, COMC 896 > ) , I WFSQ, COMC 11 16) ) • 
2(*B0W,CQM(95§34 ) , C'JH,COMC 10105 )>, C EVODD, COMC 329) ) 

1 JHK* JHCU 
XNf = NF 

DO 10 K=N1 , NF., N S 
XX * K 

XW * XK*W/XNF 
GXM = GOW#KM 
X(?IK) • XW 
GO TO 1 1 1*12 V., EVODD 

11 MF<K> • 80 M# GO SC GXW) 

GO TO 13 

12 WF<K) * BOM# 34N ( GXVO/GOW 

13 WF6QCK) = MENU **2 
XW = 2.0*H*XKyXNF 
GXH = GOHtXH 

J* K+10 

XPCJ) * XH * N 
GO TO 12,3, 4>., JHK 

2 WFCJ) * AHN# SiNHCGXH) + BHN* COSH(GXH) 

GO TO 5 

3 Wfi<J) * AHN# XH + BHN 
GO TO 5 

4 MFIJ) * AHN# SINFGXH) ♦ BHN* COStGXH) 

5 VlFSQ ( J ) * ME (*J ) **2 
10 CONTINUE 

- RETURN 
ENO 


SUBROUTINE WFX I W ( AN* , BNW, G I W, I ) 

DIMENSION XP ( 220 ) »WF( 220) , WFSQ ( 220) *JWC 100) 

COMMON COM 

EQUIVALENCE ( A,COM( 1) ) , { W,COM< 3) ) , (N.COMt 16) ) , 
l ( Nit COMC 17) ) , ( NF > COM l 10) ) f (NS, COMC 19) ) , 

2( XP,COM< 576) ) , ( ^F, COM ( 8 96) ) , ( WF SQ , COM ( 1 1 1 6 ) ) , { JW, C0M110205) ) 
1 1 = 1 

JWK = JW( I ) 

XI = 1-1 
XNF = NF 


2 

3 

4 

5 

10 


DO 10 K= N I , NF, NS 
XK = K 

XW = 2 *0*W*XK/XNF 
XPW = XI* A— W ♦ XW 
GXP W = G I w*X W 
KP = 20* ( 1-1) + K 
XP(KP) = XPW 
GO TO C2, 3, 4), JWK 
WF(KP) = ANW*SINH(GXPW) 
GO TO 5 

wF(KP) = AN W* XW + BNW 
GO TO 5 

WF(KP) = AN W* S INC GXP W ) 

WFSQlKP) = WF(KP) **2 

CONTINUE 

RETURN 

END 


+ BNW* COSH(GXPW) 


BNW* COSCGXPW) 





SUBROUTINE WFX I H ( ANH f BNH , G JH, I ) 

DI PENS IGN XP(220) ,WF( 220) ,WFSQ< 220) , JH< 100) 

COMMON COP 

COMMON COCOCM (10210) 

EQUIVALENCE (A t COM< 1) ) , (H,C0Ml(3) ) , ( H, COM ( 5 > > , ( N , COM ( 1 6 >) , 
1(NI ,COM (17) ) , ( NF *COM{ 1 8 ) J , ( NS, COM ( 19) ) , 

2 ( XP ,COM ( 576 ) } , (WF, COM (896) ),(WFSQ,COM( 1116) ), ( JH,COM( 10105) ) 

1 1=1 

JHK * JH ( I ) 

XN = N 
XI * 1-1 
XNF = NF 

CALL SLITET (2, KOOOF X ) 

GO T0(ll,6), KOOOFX 

6 DO 10 K= N I , NF , NS 
XK - K 

XH = 2 • 0 *H* XK/ X NF 
XPH = X I * A + W + XH 
GXPH = G JH *XH 
KP ^ 20*(I-1) + K +10 
XP(KP) = XPH 
20 GO TO (2,3,4), JHK 

2 WF(KP) = ANH*S I NH ( GXPH ) + BNH* COSH(GXPH) 

GO TO 5 

3 WF(KP) = ANH*XH + BNH 
GO TO 5 

4 WF(KP) - ANH* S I N ( GX PH ) + BNH* COS(GXPH) 

5 VnFSQ(KP) = WF(KP) * *2 
1C CONTINUE 

GO TO 15 

11 DO 12 K = M , NF , NS 
XK = K 

XH = 2 • 0*H* XK/ X NF 
XPH= XN * A ♦ W + XH 

GXPH = G JH * I X PH —A* ( XN + 1.0) + W) 

KP = 2 0 * ( I— 1 J + K +10 
XP(KP) = XPH 
GO TO (7,0, 9), JHK 

7 WF(KP) = ANH*SINH(GXPH) 

GO TO 13 

8 WF(KP) * ANH*XH 
GO TO 13 

9 WF(KP) = ANH* SIN(GXPH) 

13 WFSQ(KP) * WF(KP) **2 

12 CONTINUE 
15 RETURN 

END 
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SUBROUTINE PLOTFXIEV) 

C PLOTS XP VERTICALLY 

DIMENSION XPI220) ,WF(220) , WFSQ I 220) 

DIMENSION XDOWN 12 15) , YACROS { 430 ) , KKK ( 14 ) » P C 111 
COMMON COM 

EQUIVALENCE IN, COM! 16) I f ( VO, COM { /) ) , ( VOP , COM ( 325 ) ) , 
l ( XP » COM I 576) ),( WF , COM (896) ) , ( WFSQ, COM! 1 1 16 1 ) , I E VQOD, COM ( 329 ) ) f 
21 XOSIZE, COM 1807) ) , { LTM, COMt 814) ) , ( LTN,COM( 815) ) , I KN, COM (816) ) , 

3 1 K SX , CUM (817)) , I K S Y , COM { 8 1 8 ) ) , ( FX, COM! 819 ) ) , ( DX , COM ( B 20 ) ) , 

4 ( F Y , COM (821)) , ( DY, COMt 822) ) 

INTEGER XDS 126, E VODD 

1 00 10 10 = 1,5 
J = 2*10 

XDOWN I 10) ^ XP ( J ) 

YACROSI 10) = WF ( J ) 

ISO = XDSIZE +■ 10 
10 YACROSI ISQ) = WFSQ I J ) 

C WF AND WFSQ FOR ZEROTH WELL NOW STORED IN YACROS 

2 NPL = 20* I N+ 1 ) 

DO 20 II = 11, NPL 
J1 = II -5 
XDOWN I J 1 ) =XP( ID 
YACROSI Jl) = WFI ID 
JSQ = XOSIZE + Jl 
20 YACROS I JSQ) = WFSQ! ID 
C ALL WF AND WFSQ NOW STORED IN YACROS 
C KN IS THE NUMBER OF CURVES 
KKK I 1 ) = 54 

KKK ( 2 ) - KN 

C NO OF POINTS IS = TO THE VALUE OF XDSIZE 
KKK ( 3 ) = XDSIZE 
PI l) = 1.0 

C LTM SPECIFIES NUMBER LINE SPACES BETWEEN GRID LINES 
P ( 3 ) = LTM 

C LTN SPECIFIES NUMBER OF PRINT SPACES BETWEEN GRID LINES 
'■p I 4 ) = LTN 
P I 6 ) = KSX 

C FX USED TO SPECIFY STARTING VALUE OF VERTICAL SCALE 
P I 7 ) - FX 

C DX USED TO SPECIFY CHANGE IN VERTCAL GRID VALUES EACH LINE SPACE 
PIS) =0X 
PI 9) = KSY 

C FY USED TO SPECIFY STARTING VALUE OF HORIZONTAL SCALE 
PI 10) = FY 

C DY USED TO SPECIFY CHANGE IN HORIZONTAL GRID VALUES EACH PRINT SPACE 
Pill) = DY 
C TITLE 

WRITE (6,61) N,EVODO 

61 FORMAT! 2HPT,40X, 36HWAVE FUNCTION (*) AND WFSQI+) FOR N= 13, 5X, 
16HE VODD= I 2 ) 

3 CALL PLOTMYIXDOWN, YACROS, KKK, P) 

C LEGEND 

WRITE 16,62) VO , VOP , E V 

62 FORMAT I 2HPL , 30 X , 3H V0-F5. 1 , 5X, 4HVQP- 1 PEI 5. 8 , 5X , 3HE V- 1PD23 • 1 6 , 

1 3HDN= ) 

RETURN 

END 
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